Evolution of the population of Microtus Epiroticus: the 

Yoccoz-Birkeland model. 

J. J. Nieto* M. J. Pacifico* J. L. Vieitez* 
January 31, 2011 



Abstract 

We study the discretized version of a dynamical system given by a model proposed by 
Yoccoz and Birkeland to describe the evolution of the population of Microtus Epiroticus 
on Svalbard Islands, see http://zipcodezoo.eom/Animals/M/Microtuslepiroticus. We prove 
that this discretized version has an attractor A with a hyperbolic 2-periodic point p in it. 
For certain values of the parameters the system restricted to the attractor exhibits sensibility 
to initial conditions. Under certain assumptions that seems to be sustained by numerical 
simulations, the system is topologically mixing (see definition 14. 1|) explaining some of the 
high oscillations observed in Nature. Moreover, we estimate its order-2 Kolmogorov entropy 
obtaining a positive value. Finally we give numerical evidence that there is a homoclinic 
point associated with p. 
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1 Introduction 

We study the evolution of the population of Microtus Epiroticus (sibling vole) on Svalbard Islands in the Arctic 
Ocean, using a model proposed by J. C. Yoccoz and H. Birkeland, see |Ar| . It is known that there are no significant 
predation of these small mammals but in spite of that, the population presents high oscillations in its number 
albeit the lack of food is not a determinant factor to the occurrence of these phenomena. This population exhibits 
dramatic multi-annual fluctuations, by a factor greater than 20, [Ylj . 

The Sibling Vole (Microtus Epiroticus) is a species of vole found through much of northern Europe. First 
discovered in 1960 in the Grumantbyen area, they were thought to be the Common Vole until a genetic analysis 
correctly identified them in 1990, [FJASY| . 

Since these rodents were introduced from Russia on Svalbard Isles between 1930 and 1960, [YI] , the annual 
oscillations of their number may be explained, at least in part, by a non total adaptation to the environment, 
and by the pronounced seasonal fluctuation in climatic variability at Svalbard where temperatures of —30 degrees 
Celsius are common, see [YII ILBY| . 




Figure 1: Microtus Epiroticus. 



Let us first sketch the taxonomy classification of Microtus Epiroticus. 

• Domain: Eukaryota 

• Kingdom: Animalia 

• Phylum: Chordata 

• Class: Mammalia 

• Order: Rodentia 

• Family: Muridae 

• Subfamily: Arvicolinae 

• Genus: Microtus 

• Species: Microtus Epiroticus 

Jean Christophe Yoccoz and H. Birkeland, see |Arj . have proposed the following equation 

Mi 

N(t)= N(t-a)m(N(t-a))m p (t-a)S(a)da, (1) 
Ja 
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to model the evolution in time of the population of Microtus Epiroticus. In the equation it is taken into account 
only the number N(t) of fertile females at certain time t. Indeed, the inclusion in a model of the number of 
males is justified when there are difficulties for a female to find a male (for instance if the density of population 
is too small or if the ratio male-population : female-population is far away from 1 : 1) which is not the case for 
these rodents. In fact as has been pointed out by R. A. Ims in |Ims| . "spatial clumping of sexually receptive 
females induces space sharing among male voles" which implies that it is not difficult for a female to find a male. 
Moreover, the quantity of females is about the same as those of males for these rodents, |Ims21 lYl] . 

Let us describe the parameters of the model given by equation |T]): 

1. t: is the time measured in years. 

2. N(t): is the population of active females at time t. 

3. Aq: is the maturation age, 

4. Ai: is the maximal age expected for Microtus Epiroticus. 

5. m(N): annual individual reproduction rate for a population of TV individuals 

6. m p (t): is the reproduction probability at time t of the year. 

7. S(a): probability to survive up to a years. 
The model take into account the following facts: 

(a) The age when the females of Microtus Epiroticus have their first offspring is about 50 days, i.e., Aq m 0.14 
years (see [YIS| ) . 

(b) The maximal age of survival is about 2 years, i.e., A\ = 2 (see |YIj). 

(c) The seasonal factor m p (t), that is, the reproduction probability at time t of the year, varies sharply from 
in Winter to 1 from Spring to Autumn. Thus, the definition of m p (t) we adopt is 

m _ / if < t mod (1) < p 
m " [ ' ~ \ 1 if p < t mod (1) < 1 ' 

(d) The annual individual reproduction rate m(N) for a population of TV individuals, is too high when N(t) is 
small. Indeed m(TV) of the order of a constant mo > 30 individuals is realistic due to the high fertility of 
these rodents. The value of m(N(t)) decays sharply when the population N(t) increases. Following |Ar| . 
for m(TV) we adopt 

, , [ mn if TV < 1 , 
™( N ) = { moN ~, if JV"> 1 ' t >1 - 

We will assume that 7 > 1 and for some calculations we take 7 = 8.25. The reason for that is that there 
is numerical evidence, see [Ar] , that for this value of the parameter we have chaotic behavior. 

(e) Finally for the survival probability S(a), again following |Arj . we consider a linear function: 

S(a) = 1 — , if < a < Ai, and S(a) = elsewhere . 

Ai 

Remark 1.1. Another choice of functions for S(a), for instance S(a) = exp(— na), with k > 0, are also usual 
in the literature. It would be interesting to test the model given by {TJ replacing the linear function at (e) by a 
exponential one. 

Let us describe how the integral equation 

Mi 

N(t)= / N(t - a)m(N(t- a))m p (t - a)S(a)da arises. 

For N(t), the contribution of females of age in between [a, a + Aa] C [Ao, Ai] is 

fem(t — a) x (reprod. rate(£ — a) x ( season factor(t — a) x (prob. survive) x £([a,a + Aa]) 
= N(t - a) x m(N(t - a)) x m p (t - a) x S(a) x Aa , 
where £(J) is the length of the interval J and fem(t) is the number of females at time t. Here we assume that a 
female of age near Ai can reproduce and Aa is small. Taking a partition {ao = t — Ai, 01, . . . , a n =t — Ao} of 
the interval [t — Ai, t — Aq] we find 

n-l 

N(t) ~ N(t — a,j)m(N(t — aj))m p (t — aj)S(aj)Aaj , where Aa^ = (aj+i — a,j) 
3=0 

Letting n — > 00 we get at the limit the integral equation given by |T]). 



3 



1.1 The discrete model. 



There is no special reason to prefer the continuous model above to its discretization: most of the quantities 
involved, as N(t) and m(N), are by nature of discrete type. Moreover, from the experimental point of view, it 
is more natural to split the year on days and even in groups of days since it is very difficult to monitor N(t) 
experimentally. Hence, we assume that the year is split into p equal parts. 

Since the expected value of survival is bounded by A\ — 2 years we will study the evolution of N(t) for 
discrete values of t, modeling the period [0, A\] as a vector of A\p + 1 real entrances, from t — at the initial 
time of the first year, to t = 2p corresponding to Ax — 2 the final time of the second year. 

In this case the probability of survival at age £ is given by 

8{j)=l-4-, J = 0,I,2...,2p. 
2p 

where Aip = 2p. It is also convenient to consider S(j) = 1 — 2p +i • This takes into account the case where 
S(2p) > 0, i.e., when these animals can reproduce till the final of their lives. 

Given an initial vector value (N , N lt N 2 , . . . , N AlP -i, N AlP ) G M Alp+1 , the evolution of N(t) =N t ,teN, 
is governed by 

Ai p-l 

N t = N t - h m{N t - h )m p {t-h)S(h)Ah (3) 

h — A p 
2 p-l 

= - V Nt- h m(N t - h )m p {t - h)S(h) . 
p *■ — ' 

h=A p 

Next we explain the choices in equation (J3j> . 

1. We take Aop — [3^] ~ 14 which corresponds to the age at which the females have their first litter of 
pups (about 50 days). Note that if p = 100 then Ao — 0.14 corresponds to 51 days. 

2. We take Ah — - years that corresponds to the length of the unit interval in which we split the year. When 
p = 100 this gives Ah — years = 3.65 days. 

Note that the value of N at t depends only on the values of TV in \b — A\;t— Ao]. Thus, the knowledge of Nt for 
t G [— Aip, 0] (two years of observation) enables us to predict 7V t for t G [0, Aop]- When p — 100 and Ao = 0.14 
this means that the knowledge of (iVo , Ni , N2 , ■ ■ ■ , N200 ) enables us to compute /V201 , . . . , JV214. Recursively we 
may compute Nj for all j > 0. 



2 The dynamical system 

Equation Q defines a discrete dynamical system in R 2p+1 as follows: 

{N , Ni, . . . , N 2p ) h. T(N , N 1 ,...,N 2p )= (N p , N p+1 ,..., N 3p ) , 

where we have used that A\ = 2 and T : M 2p+1 — > M 2p+1 is defined recursively by equation ((3j for t = 
2p + l,...,3p. 

In order to describe theoretical properties of a system given by the discretized version Q of Yoccoz-Birkeland 
equation {l}, let us assume the following restrictions that weaken those given by conditions (a)-(e) described 
before. Doing this allows to apply the conclusions to different species respecting equation © and those restrictions. 
In particular, these conclusions will apply to the original system modeling Microtus Epiroticus. 

1. m(N) is a continuous function, m : JR + — > M + , 

2. there is mo G JR + such that 

J m > m(N) > m /2 if N < 1 and . . 

\ m N-~< > m(N) > min{^,m • A^" 7 } if N > 1 , ^ 

3. < m p (t) < 1 (so that we now allow < m p (t) < 1 for certain values of t), 

4. There is e > such that m p (t) = 1 for t in an interval of length 1 — p — e > 0, 
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5. < 2Ao < Ax and Aq + 1 < Ax (this means that in average each individual has at least two opportunities 
to reproduce), 

6. Defining Co as Co = ^ {^2,h=AoP+{p+<i)p ^C 1 )) we rec l mre c o mo > 2. From the definition of S{h) it follows 



~j( if (.-^.(.^(..(fii^^)). m 

U=A p+(p+ e )p / 



The condition Co mo > 2 will imply, as we will see below in Proposition 12.31 that the population does not 
extinguish, that is, it has the permanence property (see Definition 12, 

7. The exponent 7 satisfies 7 > 1. 

The following proposition shows that a dynamical system governed by equation §3§ and respecting the re- 
strictions 1. to 7. above is bounded. 

Proposition 2.1. For all t = 1,2, ... , Aop we have Nt < N m ax -~ m o ^ Al 2 ^°' > ^ • 

Proof. Since 7 > 1, inequalities imply Njm(Nj) < mo for all j. Moreover from m p < 1 we obtain 

A 1P -1 

N t = ^2 N{t-h)m{N(t~h))m p {t~h)S(h)Ah< 

h = A p 

A 1 p Aip A-ip 

V mo 5(ft)i = ^ V S{h) = ?± V (1-4-) = 

A P P h A P ,A P Al 

— ({Ax - A )p Alp ( AlP + ^ ~ Aop{A p + 1) 



p \ 2Axp 

Ax{Ax + l/p)-A {Ao + l/p) 



m {Ax - A ) 



2Ax 



1 1 a ,s A\-A\ \ ( {Ax - Aof 

mo {Ax - Ao) = mo 1 



2Ax J V 2Ax 



□ 



By induction we obtain that for all t > 0, N t < N max . 
Remark 2.2. For the values Ax = 2, Ao = 0.18, mo = 50 we have N max ~ 41.4 . 

2.1 Permanence. 

In this section we verify that the population given by equation Q and respecting the restrictions 1. to 7. above, 
in particular conditions (j4]) and ([5}, does not extinguish. 

Definition 2.1. We say that a system P{t) modeling the evolution of a population is permanent, or satisfies the 
permanence property, if for any positive initial vector value Po, there is e > such that the solution P{i) satisfies 

liminf P{t) > e. 
i>0 ~~ 

If a given system is permanent then, assuming that the environmental conditions do not change in time, 
the associated population will not extinguish. Thus, concerning with population dynamics this property is very 
important. 

The next proposition shows that the system under study satisfies the permanence property. 

Proposition 2.3. // i{N) = min{iV t , t g [-pAx, 0]} > then N t > for all t = 0,1, ... , A p. Moreover, 

• Ifi(N) < N^ a l then Nt > ^i{N) > i{N), t g [0,pA o ]. 

• Ifi(N) > Ni~Z then Nt > ^N^l, t g [Q,pA ]. 
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Proof. If N(t -h)<l then, by @, 

iV(t - h)m(N(t - h)) > N(t - h)™ > i(N)™ . 
Otherwise N(t — h) > 1 and then, again by Q, 

N(t - h)m(N(t-h)) > min{JV(t - fc) 1 "^, 7V(t - ft)^} > mirl {iV^^, i(7V)^} . 
Hence we have 

2p-l 

AT t = ^ N(t-h)m(N(t-h))m p (t-h)S(h)Ah> 

h=A p 

2p-l 

imin{iVi-2^,i(iV)^} £ m p (t-h)S(h)> 

, ( A p+p \ 

. r i_ 7 m ., .m \ 1 / ci ft \ ■ f Ari-7 c ° m ° ■fAr\ c o m o\ 

min|iV m J — ,t{N) — | - 2^ (! - rj") = mm|iV ma 2— — ,i(JV)— j ■ 

\h = A p+(p+e)p P 1 J 

Since, by (J5]l, Co mo > 2 we get by induction that N(t) > for all t G [0, Aop]. □ 
Clearly Proposition 12.31 implies that Nt > for all t > 0. 

Corollary 2.4. There is to > 0, depending on the initial vector value, such that we have N(t) > CQ " g N^ncH, 
t>t . 

2.2 Existence of fixed points. 

The following corollary is a straightforward consequence of Propositions 12,11 and 12.31 

Corollary 2.5. If N t > for all t G [-A p , 0] then there is t > such that ^p-N^ a Z < N t < N max for t>t . 
In particular T maps the compact set 

K = ^ c °™° N^r a Z, N max ^ into itself. 

□ 

Now set 

H:={N= (N , Nt,...,N 2p )e R 2p+l : Vj = 0, 1, . . . 2p : N 3 > 0} . (6) 

Observe that Proposition 12.31 together with Corollary 12.41 imply that T maps H into itself. 
Next we prove that T : H — > H is Lipschitz. 

Lemma 2.6. T : H — > H is a Lipschitz function. 

Proof. We put in lR 2p+1 the sup norm: ||x|| — ||(a;o,:ri, . . . ,a;2 P )|| = sup t=0 2p \ x j\- 

From the definition of T we have T(No, Ni , . . . , A^p) = (N p , N p+i , . . . , Ns p ). Hence for all j = 0, . . . ,p we 
have 

\(T(N)-T(N'M = \N j+p -N^ +p \ < \\N-N'\\. (7) 
For j — p + 1, . . . , 2p, the difference \N^ t _ h - ) m(N^ t _ h - ) ) — N l ^ t _ h ^m(N^ t _ h s ) )\ can be estimated as follows: 

(a) If Af( t _ h ) < 1 and N'i t _ h \ < 1 then by inequalities © we have that 

\N it _ h) m(N {t _ h) ) - N[ t _ h) m{N[ t _ h) )\ < m \N (t _ h) - N[ t _ h) \ . 

(b) If N( t _ h ) > 1 and > 1 then, again by Q, we have that 

\N (t . h) m(N (t . h) ) - N' (t _ h) m{N{ t _ h) )\ < \(N (t _ h) f-^ ~ (iV^o^K . 
By the Mean Value Theorem, there is TV G (TV( t _ h j, N'( t _ h \) such that 

KiV^)) 1 -- - (iV,',.,)) 1 - 7 ! = |1 - j\N-<\N^ h) - N[ t _ h) \. 
Since 7 > 1 and TV > 1 we obtain 

moKN^) 1 -^ - (N'^) 1 -^ < m ( 7 - l)|JV (t _ h) - JV ( ' t _ h) | . 
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(c) If one of the above quantities is greater than 1 and the other is not, say N[ t _ h *, > 1 and N^ t _ h - ) < 1, then 
\N it _ h) m(N {t _ h) ) - Nl t _ h) m(Nl t _ h) )\ = m \N (t _ h) - (A^.,,,) 1 - 7 ] . 
If iV(t-h) > (A^t't-h)) 1 " 7 then, since < N {t -h) < 1 and 1 - 7 < we get 

m \N (t _ h) - (N'^) 1 -^ = m (N (t _ h) - (A^,) 1 " 7 ) < 



m ((iV (t _, ) ) 1 -^(7V ( ' t _ h) ; 



l-7\ 



m \(N (t _ h) ) 1 -"< - (A^,) 1 " 7 ! < mo( 7 - l)|iV (t _ fe) - N' (t _ h) \ . 
Otherwise, if N (t _ h) < (N^^) 1 ' 1 then, since N' (t _ h) > 1 and 1 - 7 < 0, we have > N^—h) 
(A^'t-h)) 1-7 > N ( t-h) - N[ t _ h) and therefore 

\N (t - h) m(N (t _ h) ) - N' (t _ h) m(N' (t _ h) )\ = m \N (t _ h) - (N^) 1 -^ < m \N {t _ h) - N' (t _ h) \ . 

Next, to estimate \N t — AT t '| for t = p,p + 1, . . . , Aop, we use (a), (b) and (c) above as below. Let L = 
max{mo,m( 
obtain that: 



max{mo,mo(7 — 1)}. Taking into account that m p (t — h) and S(h) are between and 1 and A/i = - we 



\N t -N' t \ 



2p-l 2p-l 

N (t „ h) m{N (t „ h) )m p {t-h)S(h)Ah- ^ N{ t _ h) m{N{ t _ h) )m p {t - h)S(h)Ah 



2p~l 

- ^2 (N ( ^ k) m(N (t _ h) )-Nl t _ h) m(N' (t _ h) ))m p (t-h)S(h) 

P h=A p 
2p-l 

- \N(t- h) m{N {t _ h) )^N{ t _ h) m{N{ t _ h) )\rn p {t-h)S{h)< 



P 



h=A p 



2p-l 



2p-l 



(8) 



- L\N(t-h)-N{ t „ h) \ < - Lmax\N (t _ h) - N' {t _ h) \ < (A 1 - A )L\\N -N'\\. 

h—Aqp h—Aqp 

Taking into account that 1 < (Ai — Ao)L and inequalities ([7]) and ([8| we have that for all j = 0, . . . , p,p + 
l,...,p + A p, \T(N) - T(N')\j < (Ai - A )L ■ \\N - N'\\. By induction, since A Q p > 1, we obtain that 
\\T(N) - T(N')\\ < (A! - A Q )L ■ \\N - N'\\ finishing the proof. □ 

Corollary 2.7. There is a fixed point p for T : AC — > AC. 

Proof. By Lemma [2.6l the map T is Lipschitz hence continuous. Moreover AC is a (2p + l)-dimensional topological 
disk. Hence Brouwer Fixed Point Theorem applies, |Sp| Chapter 4, Section 7]. □ 

Remark 2.8. Since every two years (A\ — 2) the rodent population is renewed perhaps it is more natural to 
search for fixed points for T 2 : AC — > AC. So, we are interested in both, fixed points and period-two points N £ AC. 
Their existence is guaranteed by Corollary \2. 7| 

In Appendix E we estimate the coordinates of a fixed point pofT 2 :H^H. We find that the distance given 
by the norm of the supremum between p and T 2 (p) is about 8.0148 x 1CP 14 and the I 1 norm is about 4.0353 x 10 -12 . 
This estimate of p is better than that obtained by Arlot, \Ar\ Section B.8], which is of order 10 -4 for the I 1 norm. 



3 Existence of an attractor for the discrete model. 

Proposition 3.1. Let A = p| n>Q T n (AC). Then A 7^ is compact T-mvariant and there is a neighborhood 
U = (7(A) such that T(U) C U, i.e., A is an attractor for T. 

Remark 3.2. We are not assuming that A is transitive in the definition of attractor. 
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Proof. Since T(1C) C A3 we have that C„ = n™ =0 T : '(AI) is a decreasing sequence of non empty compact subsets of 
lR 2p+1 ; Co D Ci D ■ ■ ■ D C„ D ■ ■ ■ ■ Thus, by Baire Theorem, we have that A 7^ and A is compact. 
By definition of A we have 

T(A) = T(n n > T n (K.)) c n n > T n+1 (fC) c n„> T"(x:) = A, 

proving that A is T- invariant. 

Let [K] e ■= {N G M 2p+1 : dist(N,K,) < e} and e > be so small that 

[£]« C H = {N = {No, N l7 . . . , N 2p ) S R 2p+1 : Vj = 0, 1, . . .2p : AT,- > 0} . 

By Proposition 12.31 and Proposition 12.51 for all x £ [/C] e there is n(x) > such that T m ' x '(x) G AT. By continuity 
of T, see Lemma 12.61 there is U(x) a neighborhood of x contained in H such that T n ^ x \y) G JC for all 3/ 6 U(x). 
By compactness of [/C] e there is m > such that T n ([tC] t ) C A3 for all n > ni. 
Let now ?7(A) be a neighborhood of A contained in [A3] e . 

Claim 3.1. T/iere is n > smc/i i/iai T n (K.) C J7(A) /or a// n > no- 
Proof. The proof goes by contradiction. If it were not true, for all j £ IN there would exist Xj € K, and iij > n^-i, 
such that T nj (xj) £ U(A). Since A3 is compact there exists a convergent subsequence from {T' lj 1 (a; 3 -)} . g ^y. 
Without loss we may assume that {T' lj (^OljgJV itself converges to a point z £ K,. Such a point 2 cannot 
be in A since T n i {xj) U(A) for every j G IN. But, since T n+1 (A;) C T n (A3) for all n G W, we obtain 
r^ fe) G nli T h (K.). Moreover, z G nll T h QC) , otherwise there is e > such that dist(z, n^i T h (A:)) > e. 
But nlL T h (K.) D n^l + 1 r h (A:) for all j G W, so that dist(z, n^i + ! T h (A:)) > e for every I > contradicting the 
fact that T n i+ l (xj+i) — > z when I — > 00. It follows that T n {fZ) C (7(A) for all n > no, proving the claim. □ 



To conclude the proof of the proposition it is enough to verify that there is 722 > such that T n2 (U(A)) C 
[/(A). This follows from the fact that T n °(fC) C U(A) C [/(A) C [AT] E taking n 2 = no + m, thus A is an 
attractor. □ 

It is clear that the fixed point p given by Corollary 12.71 belongs to A. In [Arl Section B.8] by numerical 
methods it is found a candidate to be a fixed point. As we have pointed out above, in view of Corollary 12.71 the 
search for such a fixed point has sense. 



Remark 3.3. By Corollarv \2.5l the basin of attraction of A is the whole set H of points with positive coordinates 
(see equation (|6} ). Moreover, since K. is a disk, we can choose U(A) simply connected in the proof of Proposition 
lg.il These facts have some theoretical implications that we discuss in section [7j 

Lemma 3.4. Assume that S(2p — 1) > and m p (l) — 1. Moreover also assume that T depends smoothly on 
N = (No, A 7 !,..., N 2p ) and that 9(N Q {Nj)) ^ 0, j = 0, 1, . . . ,p. Then the differential D N T : M 2p+1 -» R 2p+1 is 
a non singular linear map. 



Proof. Let us duplicate the (p + l)-th coordinate, N p , of N = (No, ■ ■ ■ , N p , . . . , N 2p ), i.e., we write N = 
(N ,...N p ,N p ,...N 2p ) = (N (0) ,N (1) ), and consider f(N (0) ,N (1) ) = (JV (1) , AT (2) ) where N {2) = (N 2p , . . . , N 3p ). 
Thus, since the p ih -coordinate equals the (p+l) t,l -coordinate, f(N (0) ,N (1) ) is such that U p (f(N^ ), N(i)) = T(N) 
and if T is locally injective then T is locally injective too. Here lip : M 2p+2 — > M 2p+1 is the projection 

n p (a;o, . . . ,x p , Xp+i, . . . ,X2p+i) = (xo, ■ ■ ■ ,Xp-i,X p +x, ■ ■ ■ ,x 2p+1 ) . 

Taking into account that (N 2p , . . . , Ns p ) depends on (No, ■ ■ ■ , N p , . . . , N 2p ), this artifice allows us to write 
f(N (0) ,N w ) = (N (1) ,F(N m ,N (1) )), and therefore 

/ A I Id \ 

OF I BF j 



DT 

where A is a (p + 1) x (p + 1) matrix of the form 



8N lm I dNr 



A 



/ •■•0 1 \ 
•■■0 

\ ■•■ / 
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and Id is the identity (p + 1) X (p + 1) matrix. 

To prove that T is locally injective it suffice to prove that detDT 7^ 0. Hence, since det(A) = we are left 
to prove that det ( ) 7^ 0. For this we proceed as follows. Using the expression for Nt given at equation ([3} 



and denoting 



afiVjmfJVj)) 



by h(Nj) we compute g ff and find 



/ h(N )m P (l)S(2p -1) h(N 1 )m P (2)S(2p -2) 
/i(Afi)m p (l)S'(2p-l) 




h(N p )m p (p)S(p) \ 
h(N p )m p (p- l)5(p+l) 



V 





Since by hypothesis h(Nj) 7^ the thesis follows. 

Corollary 3.5. Under the hypothesis of Lemma\3.4\ we have that T : A 



ft(iV P K(l)S(2p-l) / 



A is locally injective. 



□ 



Remark 3.6. Albeit T : H H is locally injective, by Lemma \3.4\ it is not globally injective. To see this 
assume that N max > 1, 7 > 1 and that the definition of m(N) is given by equation @. // (Ao, Ai, . . . , A^p-i) = 
(Nmax, Nmax, ■ N ma x) then for t = 2p , 2p + 1, . . . , 2p + Aop we get 



2p-l 



2p-l 



1 C -v 1 . 

N t = - 22 N t -hm(Nt-h)m p (t- h)S{h) = - ^ N max m(N ma x)m p (t - h)S(h) 



h=A p 



h=A p 



2p-l 



— N, 



1-7 



m m p {t - h)S{h) . (9) 

C h=A p 

Similarly if we put (No, N\, . . . , A^p-i) = (N^J,, A m n2, • • ■ , N^2) we obtain the same values for Nt- By 
induction we get that all values are the same for t > 2p implying that T is not globally injective. 

Let us point out that: 

1. In the original model, [Arj . m(Nj) is given by equation ([2} i.e., m(Nj) = mo if Nj < 1 and m(Nj) = moN ' 
if Nj > 1. Hence N ] m(N ] ) = m Nj if Nj < 1 and Njm(Nj) = m A 1-7 if Nj > 1 implying that 



_ d{N 3 m{Nj)) _ f mo if Nj < 1 
1 dNj ~\ m (l- 7 )A r 7 7 



if A^>1 



Since 7 > 1, we have h(Nj) ^ for all Nj / 1. 

Assuming that T 2 is C 1 , Lemma r3.4| gives that the fixed point p found at Corollary 12. 7l has all its eigenvalues 
different from zero. The numerical approximation of the eigenvalues of DT p , for the estimated value of p 
obtained by [Arl Section 4.2.7] and our own estimates gives that there is a single eigenvalue of modulus 
greater than 1 which is negative, and there are A\p eigenvalues of modulus less than 1. Hence p is a 
codimension one hyperbolic fixed point of T 2 . 

The hypothesis S(2p — 1) 7^ is reasonable: otherwise one can see that for two initial vectors N = 
(Ao, JVi, • • • , A^ 2p ) and N' = (A" , Ai, • ■ • , N 2p ) with A / A" we get T(N) = T(N'). Thus the number 
Na 1p of individuals at time A\p is not affected by the first set Ao of initial individuals. In another words 
the system looses memory for a number of years less than A\ and so the actual dimension of the domain 
of T would be less than Aip + 1. 



4 Study of A for (A Q , p, 7) = (0.18, 0.30, 8.25). 

In what follows we will assume that T is smooth ( see [Arl Section 2]) and that the calculations made for the 
parameter values (0.18, 0.30, 8.25) are accurate enough to obtain that if p is the fixed point given by Corollary |2.7l 
then the eigenvalues Ai, A2, . . . , A2 P and fi of D P T satisfies |Aj | < 1 for every j = 1, . . . , 2p, and /j, ~ —3.335, in 
particular \fj,\ > Lemma 13.41 proves that p is in fact a hyperbolic fixed point with W s (p) being a codimension 

1 Arlot in [XB Section 4.2.7], obtains that fj, w -2, 29 for the parameter values (0.15, 0.30, 8.25). 
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one manifold and W u (p) an arc. Moreover, since A is an attractor, we have that W u (p) C A from which the 
fractal dimension of A is strictly greater or equal than 1. The calculations made in [Arl Section 4.2.5] give for 
this fractal dimension a value around 1.33 from which Arlot conjectures that locally the attractor is the product 
of a line by a Cantor set. 

Here we shall discuss if for the choice Ao = 0.18, p = 0.30 and 7 = 8.25 the system given by T can be 
transitiveQ 

Definition 4.1. Let f : X — > X be a continuous map defined in the topological space X . We say that the system 
defined by f is (topologically) transitive if for every pair of non-empty open subsets A, B of X there is n G 
such that f n (A)f]B 7^ 0. The dynamical system defined by f is topologically mixing if for every pair of non-empty 
open subsets A, B of X there is N > such that f n (A) n B / for all n > N . 

In |Arl Section 5] it is pointed out the interest in studying the case where the parameters are Ao = 0.18, 
p = 0.30, 7 = 8.25: it is because the numerical simulations indicates that for this parameter choice Tu is 
transitive, see [ArL Section 4.1.3, figure 12]. Moreover, in [Arl Section 4.2.7, figures 34 and 35] the geometry of 
the attractor A is depicted from the successive iterates of the local unstable manifold of the fixed point p. This 
suggests that W u (p) is dense in A. This was confirmed by the numerical simulations done by us, see figure [2] The 
next proposition shows that if the orbit of a point in W u (p) is dense in A then TU is in fact topologically mixing. 

Proposition 4.1. Let us assume that there exists xo G W u (p) such that clos(orbit + (xo)) = A that there exists a 
homoclmic point x forp that we do do not have tangencies between the stable and unstable manifold of p and that 
forward iterates by T 2 of an unstable segment s C W u (p) has diameter bounded away from zero. Then T : A — > A 
is topologically mixing. 

Proof. Observe that by hypothesis we have in particular that closW u (p) = A. Let A 7^ and B 7^ be open 
subsets of A, i.e., there are open subsets A and B of R 2p+1 such that A = .4nA and B = B(~)A. We will prove that 
there exists no such that for all n > no we have T"(A) PI B 7^ thus proving that T is topologically mixing. Since 
W u (p) is dense in A there is n 2 > such that T" 2 (x ) G A. Thus W u (p) cuts A in an arc s containing T" 2 (x ). 
Since orbit(a;o) is dense in A there exists ni > n 2 such that T ni ~ n2 (xo) G U(p) where U(p) is a neighborhood 
of p in which we may assume that we have C -linearizing coordinates, and T ni_rl2 (s) contains an arc J which 
intersects transversally Wf oc (p), this follows from the assumptions we have done. By the Inclination Lemma, see 
[PM1 Chapter 2, §7], T n (J) C 1 -approaches on compact segments of W u (p). Let v > be the radius of a ball 
contained in B. There is no > ni such that T n °(J) is i//2-dense in A and hence T n (J) is f/2-dense in A for all 
n> n . Thus T n (J) cuts B implying that T n (A) H B / for n>n . But since W u (p) C A (A is an attractor) 
we conclude that T n (A) n B 7^ for n > no proving that T is topologically mixing. 

□ 

Remark 4.2. Roughly speaking the above result means that for the parameter values Ao = 0.18, p = 0.30 and 
7 = 8.25, from the topological viewpoint we have that all possible states (No, N\ . . . , A r 2 P ) G A are visited and so a 
chaotic behavior should be expected. On the other hand, since there are fixed points like p in A if (No,Ni . . . , A^p) 
is very near p in practice we will see the same behavior for large periods of time seeming that the population of 
these rodents is in equilibria. On the other hand the hypothesis we have assumed seems to be rather strong. 

Another consequence of the density of the unstable manifold of p in A is the following (see also Remark 16. 1[) . 

Proposition 4.3. // clos(W u (p)) = A then T 2 A : A — > A is injective. 

Proof. Indeed, T 2 is injective when restricted to W u (p), for, if it were not true, there would exist x, y G W u (p) 
such that T 2 (x) = T 2 {y). But, since T 2 {p) = p it holds that W u (p) = U nelN T 2 (W?(p)) where W"(p) is the 
e-local-unstable manifold of p. Thus there is N > such that x, y G T 2N (W^(p)) and, hence, there is an arc 
7 C W u (p) with end points x and y. Applying T 2 to 7 we find a closed loop T 2 (7) contained in W u (p) which 
contradicts the fact that W u (p) is homeomorphic to IR. 

Assume now that there are x,y G A such that T 2 (x) = T 2 (y). Since T 2 is locally injective there is n > 
such that y B(x,ri) where T 2 BI ^ X \ : B(x,r\) — > H is a homeomorphism. There exists also r 2 > such that 
Tjafj/.r-o) : B(y,r2) — > H is a homeomorphism. Hence we may find V(x) C B(x,r\) a neighborhood of x and 
V(y) C B(y,r2) a neighborhood of y such that T 2 (V(x)) = T 2 (V(y)) . Since, by assumption, W u (p) is dense in 
A, there is an arc 7 C W u (p) such that has its end points x' G V(x) and y' G V(y) such that T 2 (x') = T 2 (y') 
contradicting that W u [p) is homeomorphic to IR. □ 

2 We thank Enrique Pujals for fruitful discussions on this topic. 
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We point out that the numerical simulations presented in the appendices justify that the hypothesis assumed 
in Propositions 14. 1 1 and 14.31 are reasonable. Indeed we found: 

1. If there is a homoclinic point we must have positive entropy. We estimate in Appendix A the order-2 
Kolmogorov entropy of the attractor, |Ta| . and found a positive value w 0.75. 

2. The absence of tangencies should be checked in a certain way, at least in a neighborhood of p. In algorithm 
"homclin4" presented in Appendix B, we compute the angle between the local unstable manifold Wg{p) 
and the iterate T m (£), of an arc £ C Wf{p), for m > such that T m (£) is near p, founding in all cases 
values close to ir or radians, thus W™(p) and T m (£) are almost parallel. 

3. That there is a point in W u (p) whose orbit is dense is a rather strong assumption. But when we plot the 
image of the first 1000 iterates of a single point of the local unstable manifold Wj^Jj)), projected into M 3 
we roughly recover the image of A obtained plotting all the sequences of points pseudo-randomly generated, 
see Appendix F. Moreover, in all the simulations done in algorithm "entropia3" presented in Appendix B, 
we always obtain that if N ^ N' then T 2 (N) ^ T 2 (N'), indicating that the hypothesis of the density of 
W u (p) in A assumed in Propositions 14. 1 1 and 14.31 is consistent. 

4. That forward iterates of a non trivial segment s C W u (p) have their diameters bounded away from zero also 
is rather strong. But again in all the simulations done, in particular in all runs of algorithm "homclin4", 
presented in Appendix B, we verify that this is the case. 

5. Moreover, there are theoretical results that point out that in a setting like that of this model, we cannot 
expect T 2 A to be C 1 -robustly transitive. Indeed, by construction the attractor A is contained in a simply 
connected neighborhood U C M 201 . Then by a C 1 -small perturbation we may create a sink (see |RS| for 
instance) whose basin of attraction may contain (part of) W u (p) . Nevertheless, the type of perturbations 
we can perform with T is not arbitrary and so we cannot reject a priori that for certain parameter values 
(like (Ao,p, 7) = (0.18,0.30,8.25)) the system is transitive. 

In the following subsections we check numerically the hypothesis of Propositions 14.11 and 14.31 

4.1 Estimation of the Kolmogorov Entropy of the Attractor 

As a first step to estimate the presence of chaos in A is to verify that it has sensibility with respect to initial data. 
To do so we have made computer simulations of the system given by Q with the parameter values (0.18, 0.30, 8.25). 
That A presents sensibility to initial conditions has been pointed out by Arlot, [Arl Section 4.2.6]. To test this 
property we proceed as follows: 

1. We generate M independent initial vectors jV w = (N^\N^\. . .,N^ p ), 1 < j < M. In fact what we 
have done is to generate M = 400 files with initial data chosen in a pseudo-random way. We assume that 
these 400 initial data are independent. 

2. We iterate Mimes by T 2 so that T 2l (N U) ) can be assumed, from the practical point of view, to belong to 
the attractor. The value of i that we have chosen is I = 10000 so that we are considering T 20000 (N^'). 
For simplicity of notation we still denote this iterate by N^K 

3. We add a small noise AN^ to N™> obtaining A'- 1 ' = + AN'^\ In the computer simulations we 
choose 10~ 10 < \\AN U) \\ < 10~ 8 . 

4. We specify a initial distance do and compute for every j the integer bj such that 

\\T 2l (N u) - T 2 \N U) )\\ < d , 0<i< bj, and \\T 2h ' {N ij) - T 2bj (N u) )\\ > d Q . 

We choose do =0.1 since we observe fast divergence between the orbits when this distance is achieved. 

5. In all the simulations we have done we find that bj < 80. In fact, we change the size of the perturbation 
finding that even with 10~ 18 < || A7V (j) | < 10~ 16 , the value of bj satisfies bj < 200. We conclude that 
there are numerical evidences that Tu exhibits high sensibility to initial conditions. 

As a second step to test the chaotic behavior on A we estimate its order-2 Kolmogorov entropy K giving by 
the average time for two initially near orbits of the attractor to diverge. More precisely, K is calculated from 
the average time to that is needed for two points in the attractor, which are initially within a specified maximum 
distance do, to separate until the distance between these points has become larger than do- 

The Kolmogorov entropy of an attractor can be considered as a measure for the rate of information loss along 
the attractor or as a measure for the degree of predictability of points along the attractor given an initial data. 
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In general, a positive entropy is considered as the conclusive proof that the dynamical system is chaotic. A zero 
entropy represents a constant or a regular phenomena that can be represented by a fixed point or a periodic 
attractor, [Ta| . 

Here we apply the definitions of the order-2 Kolmogorov entropy suggested by Takens in |Ta] and by Grass- 
berger and Procaccia in |GP| . see also |GP2| . According to these definitions, we will estimate the entropy from 
the average time required for two nearby distinct orbits of the attractor to diverge. 

According to Takens |Ta| and Grassberger and Procaccia |GP| . the separation of distinct nearby orbits is 
assumed to be exponential and the time interval to required for two initially nearby points to separate by a 
distance larger than do will be exponentially distributed according to 

where K is the Kolmogorov entropy, see [GP3| . For practical purposes C(to) may be transformed into a discrete 
distribution function defined as 

C(6) = e~ KhT ", with 6 = 1,2,3,..., 

where T a is the time step between two sampled data points. Given an initial pair of independent points within a 
distance do, the variable b is the number of sequential pairs of points on the attractor such that the interpoint 
distance is for the first time bigger than do. 
To estimate K we proceed as follows. 

1. We generate Z independent initial vectors — (Nq , , . . . , N^ p ), 1 < j < Z. For practical purposes 
we take for Z the same M = 400 files used to estimate sensibility to initial conditions. 

2. We iterate Mimes by T 2 so that T 2l (N U) ) can be assumed, from the practical point of view, to belong to 
the attractor. The value of I that we have chosen is I = 10000 so that we are considering t 20000 (N u) ) g A. 
For simplicity we still denote this iterate by N U) and will denote the initial N u) by T~ 20000 (N (J) ), but 
this is just a notation; we are not claiming that T is globally invertible. 

3. For each j = 1, . . . ,Z, we write in the file number j the values of 

T - 2 °°°°(iV«), N {J \ T 2 (N U) ), T 4 (N^),...,T 2044 (N W ). 

4. Given a distance d > 0, we search for pairs of vectors T h i(N u) ), T h '(N (i) ) such that \\T h i (N U) ) - 
j*i ( jyw ) || < ^ According to |STB) the value of d should be smaller than of the absolute deviation 
8n The simulations we have done give that the mean value < N > of the population N t is about 2.335 
and the average absolute deviation 



J2 \T h {N')i- < N >| w 0.97, 



400 x 2046 x 200 

thus, we take d < 0.97 x 10~ 2 (the greater value of d we have used is d = j|g). 

5. Given d, T h i {N ij) ) and T h *(iV w ) as in item [4] above, we compute the integer b = b(i, j, hi, hj) such that 

\\T 2a (T h i(N u) ))-T 2s (T hi (N (i) ))\\ < d, 0<s<6, and 

||T 26 (T^'(-'V W) ))-^ 2i '(^ fei (^ Ci) ))ll > d - 

6. Letting M = M(d) be equal to the number of distinct pairs 

T h HN {j) ), T h *{N {l) ), l<j <i< 400, 

verifying item [4] we compute 6 = S,=i bj. The program doing this task has to take care to not duplicate 
the number of times a given pair T hj (N^), T hi (N^) is computed and also to not consider as different 
strings the one starting at s = 

||T 2s (T^(AT w) ))-T 2s (T fei (Ar (i) ))l| <d, 0<s<6, and 

\\T 2b (T h i(N u) )) -T ab (T hi (N®))\\ > d, 
with the sub-strings starting at s = so > 

\\T 2a (T h * (N U) )) -T aa (T hi {N^))\\ < d, < s < s < b, and 

\\T 2b (T h *(NW))-T 2b (T h *(N®))\\ > d. 
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7. Finally we estimate the value of the entropy K of T 2 by 



K = - — In 

Ts 



1-i 
b 



where K is the maximum-likelihood estimate of the entropy K (see STB ). 

8. We repeat the items above for several values of d. Taking d m 1/100 we find more than 2000 verifying item 
[?[, while for values of d < 1/50000 the number of such pairs is too low, less than 100. More precisely, for 
d — 1/65536 = 0.0000152587890625 we find 53 strings. This is reflected in the estimate of the standard 
deviation of the entropy: for values of d too small the sample is also small and the estimation of K is less 
accurate, as one can see in Appendix A. 

To test a confidence interval for the values obtained to K we need to estimate its standard deviation. For 
this note that the standard deviation of K can be obtained from the variance of b. To do so recall, [STB] , that 

e k 

var(b) = -r—r where k — K t s . 

(e fc — l) 2 

The standard deviation in the estimate of b, computed in item 6. is given by 



a(b) = y/var{b)/M = 



M(e k - 1) 



For large values of M, a(b) will be small. In that case we can use the derivative of the function k = — ln(l — 1/6) 
in the point k = Kt 3 to estimate the standard deviation of k. 

The values obtained for the entropy of T 2 are listed in two tables in Appendix A which contain also the 
values of d we have used and those of the standard deviation ok of the entropy. For both extreme values used 
for d, namely d = 1/128 = 0.0078125 and d = 1/65536 = 0.0000152587890625 the results are less accurate, since 
0.0078125 is "too big" with respect to 0.97 « 1, and for 0.0000152587890625 there are few sample points, see 

EH. 

Nevertheless all the estimates obtained show that T 2 A has positive entropy, which implies that T|a also has 
positive order-two entropy K w 0.37. 

Thus we have strong numerical evidence that A is a chaotic attractor. 

Remark 4.4. We do not claim that we have estimated the entropy ofT^- The calculations made has to be seen 
as an indication that the model given by equation ([3]) exhibits a chaotic behavior. Rigorous proofs are needed to 
confirm our estimations. 



5 Existence of homoclinic points: numerical approach. 

In dynamical systems the presence of chaotic behavior is often associated to the existence of homoclinic points. 
We have assumed their existence in Proposition 14.11 to obtain that A is topologically mixing. Next we check 
numerically their existence. To do it we proceed as follows: 

5.1 Approximated W^ c (p). 

Due to the fact that W u (p) is one dimensional a first attempt is to try to pick a fundamental domain in Wj" c (p) 
and search by brute force if it is possible to find a candidate to be a homoclinic point there. Problem: we do not 
know precisely W^ c (p). Moreover, the value of the fixed point p is known only by an approximate value p. But 
we know that there are only one eigenvalue fi of modulus greater than 1 of DT 2 and fi is negative. Hence, since 
the other A\p eigenvalues are small in modulus, in fact all of them have modulus less than 0.5, we may assume 
that Wf oc (p) is a ^ip-dimensional disk so if we iterate p by T 2 , since fi < we have that the segment [p, T 2 (p)] 
cuts W{ oc (p) at a unique point. By the A-lemma we have that the successive iterates of [p, T 2 (p)] by T 2 converges 

Thus for numerical simulations we can take as Wj" c (p) one of these segments. In some of our simulations we 
choose [T 38 (p),T 40 (p)] as W? oc {p) and in others we take W£ c (p) as [T 30 (p),T 32 (p)]. Observe that the length of 
[T 38 (p),T 40 (p)] is less than 10" 3 and the length of [T 30 (p), T 32 [p)] is less than 10" 4 . Hence, since the mean value 
of the data is 2.335 and that the absolute deviation is 0.97 such lengths are relatively small. 
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We subdivide the chosen segments in 10000 equal parts and iterate more than 2000 times by T every point y 
of the subdivision finding the iterate T 2 -'(y) closer to p. In order to not consider misleading solutions, we discard 
the first 20 iterates and check that the orbit of y is "returning near the point p", i.e., we check that the minimum 
distance is not achieved in the 21 th iterate. Then we create a table containing the values of y and of the iterate 
of y closer to p. Not that this procedure does not prove that any of such a point y is a homoclinic point. 

5.2 Returning points. 

After this we find the value of yo and jo that minimizes dist(T 2j (y),p). In the simulations corresponding to 
WSJp) « [T 3S (p),T 40 (p)] we find that 

yo = T3S (P) + ^o( T40 (P)- T38 m and Jo =629. 

We find a suitable sub- interval Jo such that y 6 Jo C I C [T 38 (p) , T 40 (pYl , we iterate 10 times by T 2 the 
point yo and the extreme points of the segment Jo, calling them Lo and Ro □• After this we subdivide again 
T 20 (Jo) and find a small interval Ji C T 20 (Jo) around T 20 (yo) and iterate again their end points Li,Ri and 
also T 20 (yo). We continue with this procedure finding segments In C T 20 (Ih-i) and their end-points Lh,Rh till 
we arrive to the value of jo- There are cases that we cannot iterate 10 times by T 2 because distances become 
relatively large or because we cannot assume T 20 ([Lh, Rh]) to be a straight segment and in that cases we reduce 
the step size. The final step does not have to be a multiple of 10. We found that a suitable value for the length 
of the initial segment Jo is 1.122 x 10 -7 . To validate this procedure we have to check several things: 

1. control that the length of T 20 (J; 1 ) does not increase too much: we do not accept a length greater than 
10~ 4 . If the length of T 20 (Ih) is greater than 10 -4 we reduce the step used: first to 8 iterates by T 2 and 
finally by 2 iterates by T 2 . In our computations we do not need to further reduce this number of iterates. 

2. control that the segment T 20 (I h ) (or T 16 (I h ) or T 4 {I h ) in case that we have to choose a smaller step) does 
not bend too much: we require that T 20 (Ih) behaves like a straight segment. To do so we subdivide the 
segment Ik into four equal smaller segments [Lh,L' h ], [L' h , T (yo)], [T 20h (yo), R' h ], and [R' h , Rh]- Next we 
check that after 10 iterates of these intervals by T 2 , the sum of their lengths satisfies that 

T 20 ( [L h , L' h ] ) + T 20 ( [L' h , T 20h (yo )]) + T 20 ( [T 20h (y ),R' h }) + T 20 ( [R' h , R h }) 

is almost the same as the length of T 20 {[L h , R h ]). We reject any case where the quotient between both 
quantities is greater than 1.0001, reducing the number of iterates if it were necessarjQ. 

5.3 Far from tangencies. 

After computing T J0 (yo) and the corresponding points Lh and Rh for suitable ho we compute the angle 
between [Lh ,Rh ] and [T 38 (p),T 40 (p)]. We expect to have an angle close to or 180 degrees, and in fact this is 
the case in all the simulations: we obtain for the angle the value of 3.108 x 10~ 5 radians. This is an indication 
that we are not near a tangency. 

5.4 Evidence of homoclinic points. 

1. For a suitable choice of Jo = [Lo, Ro] we compute the angle between the segments [p, Lh ] and [p, Rh ]- This 
is a key point in our calculations. Before we indicate how we proceed to do so, recall that the codimension 
one submanifold Wf oc (p) of M Alp+1 locally separates Wt Alp+1 in two regions that we denote by W s ' + and 
W '-. 

On the one hand, if Lh € W a ' + and Rh € W s '~ then [Lh ,Rh ] intersects Wf oc (p) and so we have a 
homoclinic point in this segment \Lh , Rh ]- Hence, by the A-lemma the angle between successive iterates 
of the vectors [p, Lh ] and [p, Rh ] would increase up to a value close to tt. 

3 To try to subdivide the interval J C [T 38 (p),T 40 (p)] around yo of end points T 38 (p) + (5101/10000) (T 40 (p) - 
T 38 (p)) and T 38 (p) + (5103/10000) (T 40 (p) - T 38 (p)) to obtain more precision is not a good idea since forward 
iterates by T 2 of J increases their length exponentially fast. We loose any precision in the calculus after less than 
20 iterations by T 2 . 

4 In fact at the scale we have chosen this has never been the case for reasonable values of £([Lh, Rh])- 
If the number of iterates is always 10 then we get ho = [#] . 
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On the other hand, if both points are in the same region with respect to Wf oc (p), say Lh , Rh € W s ' + , 
then the segment [Lh , Rh ] will not cut Wi oc (p) and, again by the A-lemma, we have that the angle between 
successive iterates of the vectors [p, L ho ] and [p, Rh ] goes to zero when we iterate by T 2 . In this case the 
existence of a homoclinic point cannot be guaranteed. 

In the simulations we have done, see Appendix F, we obtain that for Jo of length 1.122 x 1CP 7 the initial 
angle between [p, L h(l ] and [p,R ho ] is 1.337 radians, approximately 77 degrees. For the angle between 
[p, T 2 (L ho )] and [p, T 2 (R ho )} we obtain a value of 3.011 radians which is about 173 degrees. For the angle 
between [p,T 4 (L ho )} and [p,T 4 {R ho )} we obtain a value of 3.139 radians which is about 180 degrees and 
for the angle between [p, T 6 (L ho )] and [p, T 6 (R hQ )] we obtain a value of 3.140 radians. For the subsequent 
iterates the angle diminishes slightly but up to the 14 th iterate we find that the angle is close to ty. Thus 
in that case we find evidence that a homoclinic point exists. 

2. There are choices for the length of Jo that does not lead to such evidence. Due to the exponential dilation 
in the unstable direction the behavior is rather sensible to this value. If we choose £(Io) = 1.046 x 10~ 7 , 
instead of 1.122 x 10~ 7 , we obtain at the final step that for this value both Lh and Rh belong to the 
same local connected component of lR Alp+1 \W l a oc (p). In this case we have that the angle between [p, Lh ] 
and [p,Rh ] is 1.358 x 10~ 3 radians, the angle between [p,T 2 (L ho )] and [p,T 2 (R h[l )] is 2.922 x 10~ 5 and 
the angle between [p, T 4 (L ho )] and [p, T 4 (R hQ )] is 4.582 x 10~ 6 . This indicates that both points belong to 
the same region with respect to Wi oc (p). Thus we cannot ensure the existence of homoclinic points in this 
case. 

But as we have shown above, there are choices for the length of Jo, subject to all the mentioned restrictions, 
that render numerical evidence that we in fact do have a homoclinic point associated to the fixed point p. 

In the Appendix D we give the pseudo-code of the algorithms employed to test the existence of homoclinic 
points. 

In Appendix F we show the values of the approximate homoclinic point yo £ Wjo C (p) an d the angular values 
for the iterates [p, T 2j (L ho )] and [p, T 2j ' (R h „)\ for j = 0, 1 . . . , 7. 
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6 Appendices. 

6.1 Appendix A: numerical results for the entropy. 

The following tables gives the estimation of K T i with d varying from d = 1/128 to d = 1/2048 and d varying 
from d = 1/4096 to d — 1/65536 respectively. The values of d are evenly distributed. 
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entropy estimated 


standard deviation 


d of the estimation 





71868973392930086 


0.019083966799452565 





0078125 





72369632526302973 


0.01900037121106015 





00732421875 





72434809650524738 





019137203419564205 





0068359375 





71294980785612502 





019363360739997387 





00634765625 





72339347917167310 





019254622891041776 





005859375 





72850976773398998 





019582964739043648 





00537109375 





72335976581552290 





019863263097474933 





0048828125 





72087894452782916 





020143892470160681 





00439453125 





73055620248396048 





020338205554372255 





00390625 





73804687146973547 





020719763256481047 





00341796875 





75191849797535868 





021276865360116811 





0029296875 





73746684855086799 





022038110478093715 





00244140625 





74278320582790917 





022619630761394613 





01953125 





74484051347505878 





023688490322879650 





00146484375 





75258445595582577 





025936858832676262 





0009765625 





77155883257911018 





030932730422267861 





00048828125 



entropy estimated 


standard deviation 


d of the estimation 





79719104033302477 





038944604944862542 





000244140625 





80273603232729273 





040161477427822390 





0002288818359375 





80026054504518599 





040945123242474202 





000213623046875 





80048080034198759 





042277355365906191 





0001983642578125 





78529717602382290 





043744279577575359 





00018310546875 





78083504948163185 





045414270447522060 





0001678466796875 





79617940645818047 





047304274588613750 





000152587890625 





80929969245707451 





049841776725230737 





0001373291015625 





80683413899432532 





053337136165239741 





0001220703125 





83598552255847603 





056494298338862138 





0001068115234375 





80841322520717467 





061957796245919238 





000091552734375 





86342039883772544 





068153726594861803 





0000762939453125 





85991632434641512 





076013458407264910 





00006103515625 





90006789636726776 





091761409503824907 





0000457763671875 





79158725337319783 





122667989144921260 





000030517578125 





96758402626170560 





186693997202307350 





0000152587890625 



6.2 Appendix B: description of algorithms. 

Taking into account TR , we do not care so much about the embedding dimension and use directly as vectors of 

data those given by N = (No,Ni, . . . , A^oo). 

• A first algorithm called "ratones" is used to generate 400 files named datos[i] i = 1,2, ...,400, each of 
which contains the following data: 

1. A random seed is generated to initialize a pseudo-random generator. 

2. For each i from 1 to 400 an initial vector of dimension 201 in which every component is a real number 
N h . This real number N h is in fact a floating point number of 80 bits following IEEE 754-198fJ3 
standards for the representation, calculations and manipulations of real numbers in a computer. The 
value of every element A^ for h = to h = 199 is generated calling the RANDOM function available in 
the Software Library. The value of A^oo is calculated from equation ©. N lnlt — (No, Ni, . . . , A200) 
is stored as the first value in the corresponding file datos[i]. 

3. From equation © we compute the different values of Nh for h > 201, defining in this way recursively 

6 IEEE Standard for Binary Floating-Point Arithmetic (ANSI/IEEE Std 754-1985). Also known as IEC 
60559:1989, Binary floating-point arithmetic for microprocessor systems. 
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We discard the first 9999 iterates and stored in datos[i] the following 1024 ones, 

^20000 ^j^jinit ^ ^20002 ^jyinit ^ ^22046 ^j^init \ 

• A second algorithm that we call "ratonesl" is used to perturb randomly T 20000 (N init ) in each of the 400 files 
generated by "ratones" obtaining a vector N. The random perturbations done vary from — 2 -50 « — 10 -15 
to 2" 50 w 10 -15 in each of the /i-coordinates of T 20000 (N mlt ) for h from to 199. iV 200 is computed from 
equation @. 

• The third algorithm we use, called "sensible" , computes for each i from 1 to 400 the number bi such that 
for j = to j = bi — 1 

|| T i( T 20000^mt^_ Tj -(jy)|| < 0J iiyt,, ^20000^™*^ _ yi-^jy)!! > Q J 

We use the supremum norm in the calculations since this accelerate the computations and it is clear that 
the results do not depend on the norm used. 

• Algorithm, "sensible" , also computes the mean value < b > of bi as 

j 400 

< b >= \^ bi , 

400 ^ 

i — \ 

in all the simulations done the value of bi was less than 180 and < b >~ 100. 

• The forth algorithm, "dispersion", calculates the mean value < iV' 1 ' > of data stored in the files datos[i]. 
It calculates also the mean value of all data which gives a result of < iV >~ 2.34. 

• Algorithm "dispersion" also computes the absolute average deviation 

S N = V \T h (N j )i- < N >| » 0.97. 

400 x 2046 x 200 ^ I v ; I 

• Given a value d > the algorithm "entropia3" compares the data stored in datos[j] with that stored in 
datosfi] discarding the initial vectors (only after 20000 iterates by T we assume that the vectors are in 
A). For 1 < j < i < 400 "entropia3" searches for pairs T h i(N (j) ), T hi (N {i) ) such that their distance, 
given by the norm of the supremum, is less d. "entropia3" runs 32 times generating 32 files named info[fc], 
k = 1, . . . , 32, of records each of which contains 

1. The number i of file datos[i], 

2. the number of iterates hi by T from , 

3. the value of T h ' (N (l) ), 

4. the number j of file datos[j], 

5. the number of iterates hj by T from 

6. the value of T h i (iV (j) ). 

For values of d not so small we obtain huge files infoffc], and as d decreases the size of these files decreases. 
For computational reasons we choose d max = 1/128 (corresponding to info[l] with 6, 602 KB) and d m i n — 
1/65536 (corresponding to info[32] with 196 KB). Of course the files info[fc] contain a lot of redundant 
information since if dist(T^ (iV (j) ), T hi (N (t) )) < d and also dist(T^ +i (iV (j) ), T hi+l (N (i) )) < d, with 
/ > less than the least positive value b such that dist(T h J +b (N u) ),T h * +b (N {i) )) > d, we are storing 
(j, h j: T h i (iV«); i, hi, T h '(N^)), and also (j, hj + I, T h i +l (N^); i, K + I, T hi+l (N^)). 

• Finally the algorithm "entropia4" computes the estimation of the second order entropy, K, and its standard 
deviation using the information stored in the files info[fc] and the formulas given in [STB| . 

For this we calculate for each (j, hj,T h ' (N^); i, hi,T hi (iV w )) the least positive value b such that 

dist{T h i +b {N u) ),T h * +b (N il) )) > d. 

In order not to duplicate information, once the value 6 corresponding to (j, hj ,T h i (N u) );i, hi,T hi (N (i) )) 
is calculated, we discard in this step the records (j, h'j , T h i (iV'-''); i, h[, T hi (N^)) such that hj + b> hj or 
hi + b > h'i since these should have been taken into account in the previous step. 

Remark 6.1. Although we have not taken care of the possibility that T 20000 (iV w ) = t 20000 (7V (j) ) with i / j, 
this (very rare) possibility did not occurred in any of the simulations we have done. Moreover, in accordance with 
Proposition \4.3\ in all these simulations, in particular in algorithm "entropia3", we always obtain that if N 7^ AT' 
then T 2 (N) 7^ T 2 (N'), so that the conjecture that W u (p) is dense in A is not contradicted. 



17 



6.3 Appendix C: Pseudo-code of the algorithms employed 

Here we give the pseudo code of the programs in a language close to FreePascal, the style of programming is 

procedural. 

constants used 

A0= 0.18; p = 100; Al = 2; gamma=8.25; m0=50; rho=0.30; pipa=1024; na=400; 
type of data structures used is standard, in particular "extended" means a floating point number of 10 bytes 
and "longint" or "integer" means an integer number occupying 4 bytes of memory according to the standards of 
IEEE. We also use arrays of extended or of integer and store the data in sequential files of records. 

function S(/i:integer):extended; 

{INPUT: h G Z OUTPUT: S(h) G R + } 
begin 

if (h < 0) or (h > Al * p) then S:=0 else S:= 1 - h/(Al *p+l) 

end; 

function mrho(h:integer) :extended; 

1 if < h mod 1 < p 



{ INPUT: h & Z OUTPUT: m p {h) 



elsewhere 



begin 

entrho:=trunc(rho*p); 
if ((h mod p) <entrho) then mrho:=0 else mrho:= 1; 
end; 

function eme(N:extended) :extended; 

{INPUT: N G M + OUTPUT: m(N) G R + } 
begin 
eme:=m0; lm:=N; 
if lm>l then eme:=eme*lm**(-gamma) 
end; 

procedure comienzoazar; 

begin 

randomize; semilla~maxlongint; 
end; 

procedure AZARfvar ndongint); 

{INPUT: random^eed OUTPUT: pseudo-random numberG W} 

begin 

x— random(200000); n:=x; 
end; 

function calculo(t:integer;ene:especial):extended; 

{INPUT: t 6 Z OUTPUT: N G iR 2Alp+1 } 
type 

especial = array[1..2*Al*p+l] of extended; 
begin 
lc:=0; 

for h:=floor(A0*p) to Al*p do begin 
lc:=lc+ene[t-h]*eme(ene[t-h])*mrho(t-h)*S(h ) end; 
calculo:=lc/p 
end; 

procedure eneinicial; 

{INPUT: random; OUTPUT: first vector N G R Alp+1 } 
begin 

for i:=l to Al*p do begin nhi[i]~0; rnhi[i]:=0 end; 
for i:= 1 to Al*p do begin 
AZAR(l); nhi[i]:=l +500; { we assume that at least 500 rodents are alive} 
rnhi[i]:=nhi[i]/55000 {we normalize values; N t := 1 means 55000 rodents} 
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end; 

for i:=l to Al*p do rnhaux[i]:=rnhi[i]; 
for i:=Al*p+l to 2*Al*p+l do rnhaux[i]:=0; 
rnhi[Al*p+l]:=calculo(Al*p+l,rnhaux); 
{warning: the coordinates of the vector N begin with 1 and finishes with Al*p+1} 

end; 

procedure rnhgen(ene:rentrada;var ere:rentrada); 

{INPUT: N G R Alp+1 OUTPUT: T 2 (N) G R Alp+1 } 
type 

rentrada = array[l..Al*p+l] of extended; 
begin 
t:=l;bo:=Al*p+l; 
for j:= 1 to bo do begin rnhaux[j]:=ene[j]; ere[j]:=0 end; 

for j:=bo+l to 2*Al*p+l do rnhaux[j]:=0; 
for t:=bo+l to 2*Al*p+l do begin z:=calculo(t,rnhaux); 
rnhaux[t]:=rnhaux[t]+z end; 
for i:=l to bo do ere[i]:=rnhaux[i+Al*p] 
end; 

begin {of program "ratones"} 

{INPUT: parameter values, random data 
OUTPUT: na files of data representing time series of population of Microtus Epiroticus } 

for jj:=l to na do 
begin 

rewrite(datos [jj] ) ; comienzoazar; 
writeln( 'generating datosp ,jj ,']'); 
eneinicial; rnhgen(rnhi,rnh); 
for j — 1 to 10000 do begin 
rnhv:=rnh; rnhgcn(rnhv,rnh) 
end; {20000 iterates of T: N~iT**(20000)(N)} 
for i:=l to pipa do begin 
archi [i] . numero : = ; 
for j:=l to Al*p+1 do archi[i].serie[j]:=0; 
end; 

archi[l].serie:=rnhi; archi [2] .numero: =20000; archi[2].serie:=rnh; 
for i:=3 to pipa do begin 
rnhv:=rnh; rnhgen(rnhv,rnh); {2 iterates of T each time} 
archi[i].numero:=20000+2*(i-2); archi[i].serie:=rnh; 
end; 

for i:=l to pipa do begin write(datos[jj],archi[i]); end; 
end; {of "for jj"} 
writeln('type any key to finish'); ch:= readkey; exit 
end. {of "ratones"} 



procedure AZARl(n: extended); 

{INPUT: randomjseed OUTPUT: pseudo-random numberG R} 

begin 

x:=random; n:=x-0.5; 
end; {of AZAR1} 

procedure eneperturbl; 

{INPUT: rnhi G R Alp+1 OUTPUT: rnhi + Arnhi G R Alp+1 } 

begin 

for i:=l to Al*p do begin nhi[i]:=0; end; 
for i:= 1 to Al*p do begin 
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AZARl(l); nhi[i]:=l ; 
rnhi[i]:=rnhi[i]+nhi[i]/(2**50) 
end; 

for i:=l to Al*p do rnhaux[i]:=rnhi[i]; 
for i:=Al*p+l to 2*Al*p+l do rnhaux[i]:=0; 
rnhi[Al*p+l]:=calculo(Al*p+l,rnhaux); 
end; 

begin {of program "ratonesl"} 
{INPUT: a file "datos" generated by "ratones" 
OUTPUT: a file "datosp" representing an initial small perturbation of "datos"} 

for ii: = l to na do 
begin 

rewrite(datosp[ii]); comienzoazar; 
reset(datos[ii]); xx.numero:=-l; ayuda:=true; 
while (not Eof (datos)) and (ayuda=true) do 
begin 
read(datos,xx); 

if xx.numero=0 then begin archi[l].numero:=0; archi[l].serie:=xx.serie end; 
write(xx.numero,' serie ',xx.serie[l],' — ', xx.serie[100]); 
writeln; 

if xx.mimero=20000 then begin 
rnhi:= xx. serie; ayuda:=false 
end; 
end; 

eneperturbl; rnh:=rnhi; 
for i:=2 to pipa do 

begin 
archi [i] . numero : = ; 
for j:=l to Al*p+I do archi[i].serie[j]:=0; end; 
archi[2].numero:=20000; archi[2].serie:=rnh; 
for i:=3 to pipa do 
begin 
for j:=l to 1 do 
begin 

rnhv:=rnh; rnhgen(rnhv,rnh); 
end; 

archi[i].numero: =20000+2* (i-2); archi[i].serie:=rnh; 
end; 

for i:=l to pipa do begin write(datosp,archi[i]); end; 
reset (datos); reset (datosp); 
while (not Eof (datos)) and (not Eof (datosp)) do 
begin 

read(datos,xx); read(datosp,yy); 
write(xx.numero,' serie ',xx.serie[I],' — ',yy.numero,' serie ', yy.serie[I]); 

writeln; 
end end; 

writeln ('press any key to finish'); ch: = readkey; exit 
end. {of "ratonesl"} 



function comparar (rnhx , rnhy : rentrada):longint; 

{INPUT: rnhx, rnhy € R Alp+1 OUTPUT: or 1} 
{if "comparar" =0 then ||rnft:r — rnhy\\ < tol, if 1 then > } 

begin 
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i:=l; 

cmaux:=0; {we assume that at the beginning "comparar" is 0} 
while (cmaux=0) and (i<=Al*p) do 
begin 

if (abs(rnhx[i]-rnhy[i])>=tol) then cmaux: = l; i:=i+l; 
end; 

comparar: =cmaux; 
end; 

begin {of program "entropia3"} 
{INPUT: na files generated by "ratones", 
OUTPUT: f6 files with pairs of time series di-near, i = 1, 2, ... 16; mean value < N > of N t ; 

absolute standard deviation of iV t } 
for jj:=l to na do begin 
reset(datos[jj]); j:=0; z[jj]:=0; 
while (not Eof(datos[jj])) do begin 
for i:=l to 8 do begin 
y:=0; read(datos[jj],xx); j:=j+l; 
for h:=l to Al*p do y:=y+xx.serie[h]; 
y:=y/(Al*p); z[jj]:=z[jj]+y; end 
end; {of "while not Eof"} 

z 0j] :=z [jj]/j; 

writeln('the mean value of file datos[',jj,'] is: ', z[jj]); 
end; {of " for jj"} 
prom:=0; for jj:=l to na do prom:= prom+z[jj]; 
prom:=prom/na; writeln( 'total mean value ', prom); 
writeln( 'press any key to continue'); readkey(leer); 
for jj:=l to na do begin 
reset(datos[jj]); j:=0; w[jj]:=0; 
while (not Eof(datos[jj])) do begin 
y:=0; 

read(datos[jj],xx); j:=j + l; 
for h:=l to Al*p do y:=y+abs(xx.serie[h]-z[jj]); 
y:=y/(Al*p); w[jj]:=w[jj]+y; 
end; {of "while"} 
w[jj]:=w[jj]/j; 

writeln('the absolute deviation value for datos[',jj,'] is ', w[jj]); 
end; {of "for jj" } 
dis:=0; 

for jj:=l to na do dis:=dis+w[jj]; dis:=dis/na; 
writeln('total deviation = ', dis); 
writeln('to continue press ENTER'); readln(leer); 
{we collect data} 
tol:=l/(2**(10)); 

{"tol" is what is called d in the algorithm; here we exemplify with tol « 0.001 } 

rewrite(info); 
for jj:=l to na-1 do begin 
reset(datos[jj]); 

If (not Eof(datos[jj])) then read(datos[jj],xx); {we discard the first} 
for ii:=jj+l to na do begin 
while (not Eof(datos[jj])) do begin 
read(datos[jj],xx); reset (datos[ii]); 
If (not Eof(datos[ii])) then read(datos[ii],yy); {we discard the first} 
while (not Eof(datos[ii])) do begin 
read(datos [ii],yy); 
u : =comparar (xx.serie ,yy. serie) ; 
if u=0 then {that is: \T l (N) - T\N')\\ < tol} 
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begin 

estx.numarchl:=jj; estx.numarch2:=ii; 
estx.numiterl:=xx.numero; estx.numiter2:=yy.numero; 
estx.puntol:=xx.serie; estx.punto2:=yy.serie; 
write(info,estx); 
end {of "if"} 
end {of "while not Eof(datos[ii])"} 
end {of "while not Eof(datos[jj])"} 
end {of "for ii"} end {of "for jj"} 
writeln('teclee cualquier tecla para finalizar'); ch:= readkey; exit; 
end. {of program "entropia" } 



begin {of program "entropia4"} 

{INPUT: A file with pairs of time series d;-near, 
OUTPUT: an estimation of the second order Kolmogorov-entropy K; 
an estimation of its standard deviation a^} 
rewrite(androide) ; 
for jj:=l to na do reset(datos[jj]); 
tol:=l/(2**(10)) 
begin 
base:=l; reset(info); 
while (not Eof(info)) do begin 
read(info,estx); 

if base>tope then begin writeln('error, table too small'); halt end; 
tabla[base].numarchl:=estx.numarchl; tabla[base].numarch2:=estx.numarch2; 
tabla[base].numiterl:=estx.numiterl; tabla[base].numiter2:=estx.numiter2; 
rnhl:=estx.puntol; rnh2:=estx.punto2; j:=0; 
repeat 

rnhlv:=rnhl; rnhgen(rnhlv,rnhl); rnh2v:=rnh2; rnhgen(rnh2v,rnh2); j:=j+l; 
until comparar(rnhl,rnh2)<>0; 
tabentr[base]:=j; 
if base>l then begin 
if (tabla[base].numarch2=tabla[base-l].numarch2) and 
(tabla[base-l].numiter2+tabentr[base-l]>=tabla[base].numiter2) 
then base:=base-l {overlap of data} 
else 
begin 

if (tabla[base].numarchl=tabla[base-l].numarchl) and 
(tabla[basc-l] .numiterl+tabentr[base-l] >=tabla[base] .numiterl) 
then base:=base-l {overlap of data} 
end end; 
base:=base+l; 
end; {of "while not Eof"} 

tiempos:=0; 
for i:=l to base-1 do begin 
t iempos : =tiempos+tabentr [i] ; end ; 
tiempos:=tiempos/(base-l); entropy:=-Ln(abs(l-l/tiempos)); 
writeln('the values of bj are'); 
for j:=l to base -1 do 
begin writeln('b',j,' = ',tabentr[j],' | '); end; 
writeln( 'average of bj is <b > = ',tiempos); 
writeln(' Entropy estimated is ', entropy, ', the size of the sample is ',base-l); 

rna:=base-l; 

writeln( 'standard deviation of K is :', l/(sqrt(rna)*entropy*sqrt(tiempos*(tiempos-l)))); 
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resultado— entropy; 
desvio— 1 / (sqrt (rna) *entropy*sqrt (tiempos* (tiempos- 1) ) ) ; 
writeln('to finish press any key '); ch:=readkey; 
end. {of program "entropia4"} 

6.4 Appendix D: pseudo-code of homclin4. 

Program "homclin4" 

{INPUT: a table with the candidates to be homoclinic points} 
{OUTPUT: A point in W;" c (p) such that near it there is numerical evidence that it exists a homoclinic point} 
This program uses, apart from the functions and procedures defined above, two functions "distl2" and " angulo" . 
"distl2" computes the Euclidean distance between points, and "angulo" computes the angle between a pair of 
vectors, "angulo" uses a function "prodint" that calculates the inner product between vectors. The program 
also uses two procedures, "minimo" that computes the minimum between real data stored in a file called 
"candihomclin" and "iterar" that iterates the function T 2 a prescribed number of times. 

function distl2(rnhx,rnhy: rentrada):extended; 

{calculates euclidean distance between points} 
var cmaux,i:longint; raux,dist:extended; maximo:extended; rnhd:rentrada; 

begin 

i:—2; maximo:=abs(rnhx[l]-rnhy[l]); 
while i<=Al*p+l do begin 
if abs(rnhx[i]-rnhy[i])> maximo then maximo:=abs(rnhx[i]-rnhy[i]); 
i:=i+l end; 
if maximoOO then 
for i:=l to Al*p+1 do rnhd[i]:= abs(rnhx[i]-rnhy[i])/maximo; 
i:=l; dist:=0; { assume distance is 0} 
if maximoOO then 
while (i<=Al*p+l) do begin 
dist:=dist+rnhd[i]*rnhd[i]; i:=i+l; end; 
distl2:=maximo*sqrt(dist); 
end; 

function prodint(rnhx,rnhy: rentrada):extended; 

{computes inner product of vectors} 
var idongint; prod:extended; rnhd,rnhe:rentrada;maximox,maximoy:extended; 

begin 

i:=2; maximox:=abs(rnhx[l]);maximoy:=abs(rnhy[l]); 

while i<=Al*p+l do begin 
if abs(rnhx[i]) > maximox then maximox:=abs(rnhx[i]); 
if abs(rnhy[i]) > maximoy then maximoy:=abs(rnhy[i]); 
i:=i+l end; 
if (maximox*maximoyOO) then begin 
for i:=l to Al*p+1 do 
begin rnhd[i]:= rnhx[i]/maximox; rnhe[i]:=rnhy[i]/maximoy end; 
i:=l; prod:=0; 
while i<=Al*p+l do begin 
prod:=prod+rnhd[i]*rnhe[i]; i:=i+l end; 
pro dint : = prod * maximox *maximoy ; end 
else prodint:=0; 
end; 

function angulo(rnhx,rnhy:rentrada):extended; 

var equis, ye, zeta:extended; 
begin 

zeta:= prodint (rnhx,rnhy); equis:=sqrt(prodint(rnhx,rnhx)); ye:=sqrt(prodint(rnhy,rnhy)); 
if (equis=0) or (ye=0) then angulo:=0 else angulo— arccos(zeta/(equis*ye)); 

end; 
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procedure minimo; 

var minaux:extended;fijmin:rentrada; seguir:boolean; 
begin 

min:=l; seguir:=true; {"min" is set to a value which will not be the minimum} 
while (not Eof(refcandihomclin)) and (seguir=true) do begin 
read(refcandihomclin,homocl); if homocl.punto[2]<>0 then begin 
seguir:=false; min: = 10**(-2); refhomocl:=homocl end; 
if (homocl.punto[2]=0) and (homocl.punto[l]<min) and (homocl.numh<700) 
and (homocl.numh>10) then begin 
min:=homocl.punto[l]; refhomocl:=homocl end; 
end; {of while} 
end; 

procedure iterar(paso:integer;sota:rentrada; var sotal:rentrada); 

var rnhj,rnhjv:rentrada; {"paso" controls the number of iterations} 
begin rnhjv:=sota; 
for j: = l to paso do 
begin rnhgcn(rnhjv,rnhj); rnhjv:=rnhj; end; 
sotal:=rnhj; 
end; 

begin {of homclin4} 
while not Eof(candihomclin) do begin 
read(candihomclin,homocl) ; 
if (homocl.punto[2]=0.0) and (homocl.punto[l]j0.001) then 

write(refcandihomclin,homocl) ; 
if homocl.punto[2]j^0.0 then write(refcandihomclin,homocl) 
end; {of while} 
reset (refcandihomclin) ; 
minimo; 

writeln('minimum distance to p is ', min); 
writeln( 'value of i=',refhomocl.numi,' iterate closest to p is ',refhomocl.numh+10); 
writeln( 'initial approximation to candidate to homoclinic point M is '); 
for j:=l to Al*p+1 do begin 
fijol2[j]:= fijo6[j] + (10000-refhomocLnumi)*fijo8[j]; 
if (j mod 3=0) then writeln(fijol2[j],'— ') 
else write(fijol2[j],' — '); end; 
tolerancia: {a label of reference} 
if (rcfhomocl.numh mod 2 = 0) then techo:=refhomocl.numh+10 
else techo:=refhomocl.numh+10; 
fijo8:=restar(fijo6,fijo4); 
for j:=l to Al*p+1 do fijo8[j]:=fijo8[j]/10000; 
writeln; 

writeln('Next we refine the choice, in particular we find L and R'); 

writcln('points in [T 38 (p), T 40 (p)] identified with W£ c (p)'); 
writeln('such that M is between them and such that the iterates'); 
writeln('of L and R are in different components with respect to'); 
writeln('the local stable manifold Wi oc (p) of p.'); 
writeln('For convenience we continue to denote by M, L and R their iterates by T 2 '); 
writeln('Enter gap distance as a real exponent of 2 not greater than 30'); 
writeln('the gap distance will be 2^~ exponmt) '); 
write('To finish the program enter exponent=0, exponente = '); 
readln(scmillon); 
if semillon< then 
begin semillon:=-semillon; 
writeln('a negative value has been entered, ',semillon,' will be assumed'); 

end; 
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if semillon> 20 then 
begin 

writeln( 'exponent too large, a value of 10 will be assumed'); 
semillon:=10 end; 
while semillonO do 
begin 

tol:=2**(semillon) ; writeln('tol= ' , 1/tol) ; 
for j:=l to Al*p+1 do 
begin 

fijol2[j]:= fijo6[j] + (10000-refhomocl.numi)*fijo8[j]; 
fijoll[j]:=fijol2[j]-fijo8[j]/tol; 
fijol3[j]:= fijol2[j]+fijo8[j]/tol; 
fijoll5[j]:=fijol2[j]-fijo8[j]/(2*tol); 
fijol35[j]:=fijol2[j]+fijo8[j]/(2*tol); 
end; 

while techo> do 
begin 

if techo> = 10 then begin 
iterar(10,fijol2,rnhl); iterar(10,fijoll,rnh0); 
iterar(10,fijol3,rnh2); iterar(10,fijoll5,rnh05); 
iterar(10,fijol35,rnh25); techo:=techo-10; 
writeln('distance between left and right iterates L and R is ',distl2(rnh0,rnh2)); 

writeln('dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= ', 
distl2(rnh0,rnh05)+distl2(rnh05,rnhl)+distl2(rnhl,rnh25)+distl2(rnh25,rnh2)); 

if (distl2(rnh0,rnh2) > 0.0001) or 
(distl2(rnh0,rnh05)+distl2(rnh05,rnhl)+distl2(rnhl,rnh25)+distl2(rnh25,rnh2) 

> 1.001*distl2(rnh0,rnh2)) 

then begin 

writeln( 'distance between iterates is too large or curvature is big'); 
techo:=techo+10; iterar(8,fijol2,rnhl); 
iterar(8,fijoll,rnh0); iterar(8,fijol3,rnh2); 
iterar(8,fijoll5,rnh05); iterar(8,fijol35,rnh25); 
techo:=techo-8; 

writeln('iterating 8 times the new distance between L and R is ',distl2(rnh0,rnh2)); 
writeln('dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= ', 
distl2(rnh0,rnh05)+distl2(rnh05,rnhl)+distl2(rnhl,rnh25)+distl2(rnh25,rnh2)); 

if (distl2(rnh0,rnh2) > 0.0001) or 
(distl2(rnh0,rnh05)+distl2(rnh05,rnhl)+distl2(rnhl,rnh25)+distl2(rnh25,rnh2) 

> 1.001*distl2(rnh0,rnh2)) 

then begin 

writeln('distance between iterates continues to be too large or curvature is big'); 
techo:=techo+8; iterar(2,fijol2,rnhl); 
iterar(2,fijoll,rnh0); iterar(2,fijol3,rnh2); 
iterar(2,fijoll5,rnh05); iterar(2,fijol35,rnh25); 
techo:=techo-2; 

writeln( 'iterating 2 times the new distance between L and R is ',distl2(rnh0,rnh2)); 
writeln('dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= ', 
distl2 (rnh0,rnh05) +distl2 (rnh05,rnhl ) +distl2(rnhl ,rnh25) +distl2 (rnh25,rnh2) ) ; 

end {of inner "if then" } 
end end {of outer "if then" } 
else begin { now "techo" is less or equal than 10} 
iterar(techo,fijol2,rnhl); iterar(techo,fijoll,rnhO); 
iterar(techo,fljol3,rnh2); techo:=0 
end; 

fijo8:=restar(rnh2,rnh0); fijol2:=rnhl; 
for j:=l to Al*p+1 do 
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begin 

fijoll[j]:=fijol2[j]-fijo8[j]/tol; fijol3[j]:=fijol2[j]+fijo8[j]/tol; 
fijoll5[j]:=fijol2[j]-fijo8[j]/(2*tol); fijol35[j]:=fijol2[j]+fijo8[j]/(2*tol); 

end; 

writeln(' iterates^ ',refliomoclmumli-|-10-techo,' 
distance between endpoints L and R previous to iteration is '); 
WTiteln(distl2(fijoll,fijol3)); 
if distl2(fijoll,fijol3) i 0.00001 then for j:=l to Al*p+1 do 

begin 

fijoll[j]:=fijol2[j]-fijo8B]/(16*tol); fijol3[j]:=fijol2[j]+fijo8[i]/(16*tol) 

end; 

if distl2(fijoll,fijol3)i=0.0000000000000001 then begin 
writeln('tol is too small, please, reduce the exponent'); goto tolerancia; end; 

end; {del while techo} 
dist0:=distl2(fijo,rnh0); distl:=distl2(fijo,rnhl); 
dist2:=distl2(fijo,rnh2); writeln; 
writeln(' Euclidean dist from p to original point T 2 ',*(refhomocl.numh+10),'is ',distl); 
writeln(' Euclidean dist from p to left point ',dist0); 
writeln(' Euclidean dist from p to right point ',dist2); 
writeln(' Euclidean dist between left and right points is '); writeln(distl2(rnh0,rnh2)); 
if (rcfhomocl.numh mod 2 = 0) then techo:=refhomocl.numh+10 
else techo:=refhomocl.numh+10; 
fijol4:=restar(fijo4,fijo2); fijol6:=restar(rnh2,rnh0); rnhgcn(fijo6,fijo8); 
writeln('angle between W"(p) and iterated arc LM is = '); 
write(angulo(fijol4,fijol6)); 
writeln(' angle in degrees is approx = ' ,round(angulo(fijol4,fijol6)*180/Pi)); 
writeln(' Euclidean dist between left end-point of W"(p) and L is '); 

writeln(distl2 (fijo6,rnh0) ) ; 
writeln(' Euclidean dist between left end-point of W"(p) and R is '); 

writeln( distl2(fijo6,rnh2)); 
writeln(' Euclidean dist between right end-point of W"(p) and L is '); 

writcln(distl2 (fljo8,rnh0)) ; 
writeln(' Euclidean dist between right end-point of W"(p) and R is '); 
writeln (distl2 (fijo8 ,rnh2) ) ; 
rnhgen(rnh0,rnh0v) ;rnhgen(rnh2,rnh2v) ; 
writeln('rate of dist between rnhO, rnh2 and their iterates by T 2 is '); 
writeln(distl2(rnh0v,rnh2v)/distl2(rnh0,rnh2)); 
fijol8:=restar(fijo,rnh0) ; fijo20:=restar(fijo,rnh2) ; 
writeln('angle between vectors (p,L) and (p,R) is '); 
write(angulo(fijol8,fijo20)); 
writeln(' angle in degrees is approx = ',round(angulo(fijol8,fijo20)*180/Pi)); 
fijol8:=restar(fijo,rnh0v); fijo20:=restar(fijo,rnh2v); 
writeln('angle between vectors (p,T 2 (L)) and (p,T 2 (7?)) is '); 
writ c ( angulo (fij o 1 8 , fij o20 ) ) ; 
writeln(' angle in degrees is approx = V oun d(angulo(fijol8,fijo20)*180/Pi)); 

for ii:=l to 6 do begin 
newnx[2*ii-l]:=rnh0v;newfix[2*ii]:=rnh2v; 
rnhgen(newfix[2*ii-l] ,rnh0v) ;rnhgen(newfix[2*ii] ,rnh2v) ; 
fijol8:=restar(fijo,rnh0v); fijo20:=restar(fijo,rnh2v); 
writeln('angle between vectors (p,T 2 ',*(ii+l),'(L)) and (p,T 2 ',*(ii+l),'(R)) is '); 
writ c ( angulo (fij o 1 8 , fij o20 )) ; 
writeln(' angle in degrees is approx = V oun d(angulo(fijol8,fijo20)*180/Pi)); 

end; writeln; 

nuevofijo2:=restar(nuevofijo,rnh0); nuevofijo4:=restar(nuevofijo,rnh2); 
writcln('angle between vectors (p,L) and (p,R) is '); 
write(angulo(nuevofijo2,nuevofijo4)); 
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writeln(' angle in degrees is approx = ',round(angulo(nuevonjo2,nuevofijo4)*180/Pi)); 
writeln; fijo8:=restar(fijo6,fijo4); 
for j— 1 to Al*p+1 do fijo8[j]:=fijo8[j]/10000; 
writeln; writeln( 'Enter gap distance as a real exponent of 2 no greater than 30'); 
writeln('last exponent used is ',semillon); 
write('To finish the program enter exponent=0, exponent = '); 
readln(semillon) ; 
if semillon< then begin 
semillon: =-semillon ; 
writeln('a negative value has been entered, ', semillon,' will be assumed'); 

end; 

if semillon> 20 then begin 
writeln('exponent too large, a value of 10 will be assumed'); 
semillon: =10; 
end end; 
write(candihomclin2,homocl); 
writeln('To continue press ENTER'); read(leer); writeln( 'Press any key to finish'); ch:= readkey; 

end. 

6.5 Appendix E: coordinates of fixed point. 

Approximate coordinates of the fixed point p € Ft 201 of T 2 (N) are given in the following table. 



1.2326490487970465E + 0000 
1.171 7713776064593E + 0000 
1.1262379872047533E + 0000 
1.1049811591622351E + 0000 
1.1220867468556205E + 0000 
1.1986450973703732E + 0000 
1.1754839367637863E + 0000 
1.0755926342084036S + 0000 
9.7570133165302088S - 0001 
8.7581002909763817£ - 0001 
7.7591872654225546S - 0001 
6. 87778965 12033852S - 0001 
6.1568060751297437S - 0001 
5.5789143427761938S - 0001 
5.1301221031245340S - 0001 
4.7994290586969646£ - 0001 
1.1774339903181779S + 0000 
2.1271151896510347S + 0000 
2.9569342453691580S + 0000 
3.6876953663770639S + 0000 
4.3389071105434127S + 0000 
4.9281739025629719S + 0000 
5.5692297788707215S + 0000 
5.5004355530098589S + 0000 
5.4018415509009528S + 0000 
5.3021425008934430S + 0000 
5.2023019446101286S + 0000 
5.1024293056141891S + 0000 
5.0025458322320211S + 0000 
4.9026609256467161S + 0000 
4.8027768933916233S + 0000 
4.7028939124217354S + 0000 
4.6030121628245493S + 0000 
4.5031318571 107021S + 0000 
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1.1906886685005045S + 0000 
1.1392463406313234.E + 0000 
1.1086283238645345£ + 0000 
1.1109795220165019.E + 0000 
1.1649033774342367.E + 0000 
1.2087810376155805S + 0000 
1.1088897350601978S + 0000 
1.0089984325048151£ + 0000 
9.0910712994943241.E - 0001 
8.0921582739404970S - 0001 
7.1528062031461468S - 0001 
6.3804515773116380£ - 0001 
5.7565740489494235S - 0001 
5.2661124986901219S - 0001 
4.8971023744324276S - 0001 
8.3241286907550863S - 0001 
1.8248209517576603E + 0000 
2.6923595086534992£ + 0000 
3.4538676455028848S + 0000 
4.1295537791098290S + 0000 
4.7376142351332896S + 0000 
5.4622938956939894£ + 0000 
5.5318006317643429S + 0000 
5.4349605692581039S + 0000 
5.3354035019177187S + 0000 
5.2355878340449946S + 0000 
5.1357218479879554S + 0000 
5.0358409167684584S + 0000 
4.9359558081745136S + 0000 
4.8360714621823154S + 0000 
4.7361881120565536E' + 0000 
4.6363059301899531S + 0000 
4.5364251171897114S + 0000 
4.4365458667634358S + 0000 
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From the analytic expression of T, it is clear that T(p) 7^ p so that p has period 2. 



6.6 Appendix F: homoclinic points search. 

We plot a projection of the attractor in three-dimensional space (averaging some coordinates at the beginning of 
the year, in the middle of the year and in Spring). For that purpose we use MATLAB© . When we plot the 
image of the first 1000 iterates of a single point of the local unstable manifold W ; " c (p) we roughly recover the 
image of A obtained plotting all the sequences of points pseudo-randomly generated. This is an indication that 
W u (p) may be dense in A. The small red circle in the figures, indicates the approximate position of the fixed 

point p. 

We also give approximate coordinates of the homoclinic point yo and angular values obtained with an initial 
length of Jo of 1.122 x 10~ 7 in the appendix below with two copies of runnings of "homclin7.exe" which is a 
refinement of "homclin4.exe" which generates an output close to LaTeX. In these runs we use three values for 
the parameter "exponent", one of them is 10.80 and the other is 15.03. We also exhibit one exponent, 11.00 
which fails to detect homoclinic points. Observe that 10.80 is not very far apart from 11.00. 
We only exhibit samples of the runs since they are rather extensive. It is possible to observe that the program 
corrects the quantity of iterations when the results are larger than certain bounds. 

Runs of "homoclin7" 

Enter gap as an exponent of 2 not greater than 20 and greater than 3, we choose gap=2 exponent . 
This gap will be used to divide the distance between three consecutive points 
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Figure 2: A thousand iterates of y € W^ c {p) projected in 3D (in black) versus the projection of 
A (in green). 

of the initial subdivision of [T 38 (p), T 40 (p)] centered around 

the rough homoclinic point yo previously found. 
To finish the program enter exponent=0, 
exponent = 1.0800000000000000E+0001 gap = 1782.8875536 exponent = 10.8000000 
iter= 0, dist(L,R) previous to iteration is 3.5890677421214318E-0010 
dist between L and R after applying T 20 is 6.3294602762248157E-0007 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 6.3294602764511862E-0007 
iter= 10, dist(L,R) previous to iteration is 7.1002349673070908E-0010 
dist between L and R after applying T 20 is 1.0395709474479636E-0008 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 1.0395709474479637E-0008 
iter= 20, dist(L,R) previous to iteration is 1.1661654789708924E-0011 
dist between L and R after applying T 20 is 6.7123467133053106E-0010 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 6.7123467133053267E-0010 



iter= 600, dist(L,R) previous to iteration is 1.7874017013569691E-0011 
dist between L and R after applying T 20 is 3.5203152021803773E-0006 

dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 3.5203152021804383E-0006 
iter= 610, dist(L,R) previous to iteration is 3.9490041813784173E-0009 
dist between L and R after applying T 20 is 4.3191750633612666E-0005 

dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 4.3191750633619148E-0005 
iter= 620, dist(L,R) previous to iteration is 4.8451457912339789E-0008 
dist between L and R after applying T 18 is 2.7846424585034373E-0003 

dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 2.7846424590282723E-0003 
iter= 629, dist(L,R) previous to iteration is 3.1237443470096728E-0006 

Sup distance from fixed point p to point T 1258 (y ) is 0.00022627982 
LI distance from fixed point p to point T 1258 (j/ ) is 0.01670265394 
Euclidean distance from fixed point p to point T 1258 (y ) is 0.00141044319 
Euclidean distance from fixed point p to point L is 0.00007862956 
Euclidean distance from fixed point p to point R is 0.00280177267 
Euclidean distance between L and R is 0.00278464246 
angle between W^ip) and iterated arc LR = 0.00003 radians, angle in degrees is ~ 
angle between vectors (p,L) and (p,R) is 1.33746 angle in degrees is w 77 
angle between vectors (p,T 2 (L)) and (p,T 2 (R)) is 3.01126 radians, angle in degrees is w 173 
angle between vectors (p,T 4 (L)) and (p,T 4 (R)) is 3.13900 radians angle in degrees is « 180 
angle between vectors (p,T 6 (L)) and (p,T 6 (R)) is 3.14031 radians angle in degrees is « 180 
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Coordinates of initial approxination to hotioclinic point $y_0f ara 



1,232652605 
1,139246960 
1,105533061 
1,198633030 
1,108875395 
0,942389542 
0,775903690 
0,638029941 
0,541542333 
0,479928177 
1,824779645 
3,210446691 
4,338801024 
5,462154712 
5,467800406 
5,302039560 
5,135622250 
4,969154518 
4,802683963 
4,636216333 
4,469752522 
4,303292980 
4,136838552 
3,970390963 
3,803916723 
3,637434550 
3,470952377 
3,304474833 
3,138007377 
2,971554555 
2,805350028 
2,639336191 
2,473540491 
2,308086547 
2,143181680 
1,979185450 
1,816743387 
1,657059138 
1,502470723 
1,357707488 
1,232637195 



1.211051279 
1.126237833 
1.110974047 
1.242064077 
1.075578224 
0.909092372 
0.744619505 
0.615665389 
0.526596191 
0.471379141 
2.127067226 
3.453786886 
4.541277853 
5.569115435 
5.434854900 
5.268766684 
5.102330375 
4.935860211 
4.769390177 
4.602923232 
4.436460270 
4.270001638 
4.103548422 
3.937102462 
3.770620289 
3.604138116 
3.437656223 
3.271180416 
3.104715460 
2.938266307 
2.772133937 
2.606156503 
2.440417232 
2.275052928 
2.110294376 
1.946544654 
1.784532649 
1.625628368 
1.472515186 
1.330661505 



1.190691205 ! 
1.115883924 ! 
1.122079841 I 
1.208766906 ! 
1.042281054 ! 
0.875795201 ! 
0.715265478 ! 
0.594888025 ! 
0.512997217 ! 
0.832391559 ! 
2.415966274 ! 
3.687608184 ! 
4.737495899 ! 
5.559124568 ! 
5.401736578 ! 
5.235486232 ! 
5.069037476 ! 
4.902565995 ! 
4.736096515 ! 
4.569630292 ! 
4.403168180 ! 
4.236710511 I 
4.070258579 ! 
3.903806027 ! 
3.737323854 ! 
3.570841681 ! 
3.404360368 ! 
3.237886428 ! 
3.071424167 ! 
2.905030261 ! 
2.738923947 ! 
2.572986176 ! 
2.407308686 ! 
2.242043015 ! 
2.077446445 ! 
1.913971450 ! 
1.752442321 ! 
1.594421010 I 
1.442992241 ! 
1.304486305 ! 



1.171773334 
1.108626365 
1.139728196 
1.175469736 
1.008983883 
0.842498031 
0.687763786 
0.575642216 
O.5007O8233 
1.177406005 
2.692298311 
3.912606440 
4.928049563 
5.531692464 
5.368538383 
5.202201011 
5.035743319 
4.869271877 
4.702802982 
4.536337520 
4.369376262 
4.203419612 
4.036969042 
3.870509592 
3.704027419 
3.537545246 
3.371064836 
3.204592901 
3.038133548 
2.871798613 
2.705720592 
2.539826078 
2.374216295 
2.209059275 
2.044642259 
1.881473880 
1.720487792 
1.563467758 
1.413965467 
1.279316068 



1,154533228 
1,104978147 
1,164893175 
1,142172565 
0,975686713 
0,809200860 
0,662041411 
0,557876277 
0,489695409 
1,508319728 
2,956866483 
4,129453928 
5,113517681 
5,500328368 
5,335299888 
5,168912678 
5,002448902 
4,835977865 
4,669509585 
4,503044927 
4,336584525 
4,170128954 
4,003679830 
3,837213158 
3,670730985 
3,504248812 
3,337769649 
3,171299871 
3,004843657 
2,838571735 
2,672524458 
2,506677163 
2,341141663 
2,176104467 
2,011886743 
1,849061068 
1,688686670 
1,532804025 
1,385508654 
1,255306356 



Euclidean distance fran Sp$ to $a_0$ = 0.000894332 



Figure 3: Coordinates of (rough) homoclinic point yo associated to p. 



angle between vectors (p, T 8 (L)) and (p,T 8 (R)) is 3.13743 radians angle in degrees is « 180 
angle between vectors (p, T 10 (L)) and (p,T 10 (R)) is 3.12545 radians angle in degrees is ~ 179 
angle between vectors (p, T (L)) and (p, T 1 (R)) is 3.10750 radians angle in degrees is £s 178 
angle between vectors (p,T 14 (L)) and (p, T 14 (R)) is 3.05523 radians angle in degrees is £s 175 

Enter gap distance as a real exponent of 2 between 3 and 20 
last exponent used is 10.80000000 
To finish the program enter exponent=0. Chosen exponent = 15.03200000 
gap = 33502.9380910 exponent = 15.0320000 
iter= 0, dist(L,R) previous to iteration is 1.9099531465388560E-0011 
dist between L and R after applying T 30 is 1.4460501331114946E-0008 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 1.4460501331114946E-0008 
iter= 15, dist(L,R) previous to iteration is 8.6323753141429872E-0013 
dist between L and R after applying T 30 is 1.6945248777241631E-0009 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 1.6945248777241637E-0009 
iter= 30, dist(L,R) previous to iteration is 1.0115684115597356E-0013 
dist between L and R after applying T 30 is 2.4222230946916936E-0009 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 2.4222230946916942E-0009 
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iter= 594, dist(L,R) previous to iteration is 1.2923688152099415E-0009 
dist between L and R after applying T 30 is 9.5328543032789122E-0004 

dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 9.5328543041795548E-0004 
distance between iterates is too large or curvature is big 
dist between L and R after applying T 16 is 6.3306020056154999E-0008 

dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 6.3306020056154999E-0008 
iter= 602, dist(L,R) previous to iteration is 3.7791322627648531E-0012 
dist between L and R after applying T 30 is 1.9472974658090783E-0005 

dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 1.9472974658261263E-0005 
iter= 617, dist(L,R) previous to iteration is 1.1624636981219546E-0009 
dist between L and R after applying T 24 ) is 2.7688866606693677E-0003 

dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 2.7688866611853443E-0003 
iter= 629, dist(L,R) previous to iteration is 1.6529216948955862E-0007 

Sup distance from fixed point p to point T 1258 (j/ ) is 0.00022627982 
LI distance from fixed point p to point T 1258 (y ) is 0.01670265394 
Euclidean distance from fixed point p to point T 1258 (j/ ) is 0.00141044319 
Euclidean distance from fixed point p to point L is 0.00279389720 
Euclidean distance from fixed point p to point R is 0.00008060930 
Euclidean distance between L and R is 0.00276888666 
angle between W"(p) and iterated arc LR — 3.14156 radians, angle in degrees is « 180 
angle between vectors (p, L) and (p, R) is 1.24158 angle in degrees is ^ 71 
angle between vectors (p,T 2 (L)) and (p,T 2 (R)) is 2.76672 radians, angle in degrees is ~ 159 
angle between vectors (p, T 4 (L)) and (p,T 4 (R)) is 3.13482 radians angle in degrees is ~ 180 
angle between vectors (p, T 6 (L)) and (p,T 6 (R)) is 3.14036 radians angle in degrees is ~ 180 
angle between vectors (p,T 8 (L)) and (p,T 8 (R)) is 3.13746 radians angle in degrees is « 180 
angle between vectors (p, T 10 (L)) and (p, T 10 (R)) is 3.12554 radians angle in degrees is £s 179 
angle between vectors (p,T 12 (L)) and (p,T 12 (R)) is 3.10770 radians angle in degrees is ~ 178 
angle between vectors (p, T 14 (L)) and (p, T 14 (R)) is 3.05600 radians angle in degrees is ^ 175 

Enter gap distance as a real exponent of 2 between 3 and 20 
last exponent used is 15.03200000 
To finish the program enter exponent=0. Chosen exponent = 11.00000000 
gap = 2048.0000000 exponent = 11.0000000 
iter= dist(L,R) previous to iteration is 3.1244649405582926E-0010 
dist between L and R after applying T 24 is 7.0403463016490730E-0008 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 7.0403463017652905E-0008 
iter= 12, dist(L,R) previous to iteration is 6.8753381889103818E-0011 
dist between L and R after applying T 24 is 1.6655785942826604E-0007 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 1.6655785942826604E-0007 
iter= 24, dist(L,R) previous to iteration is 1.6265415954054225E-0010 
dist between L and R after applying T 24 is 4. 167433222936225 1E-0008 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 4.1674332229362251E-0008 
iter= 36, dist(L,R) previous to iteration is 4.0697590197929017E-0011 
dist between L and R after applying T 24 is 3.5735658309920103E-0007 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 3.5735658309920103E-0007 
iter= 48, dist(L,R) previous to iteration is 3.4898103807564854E-0010 
dist between L and R after applying T 24 is 3.1598931331358240E-0006 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 3.1598931331363929E-0006 
iter= 60, dist(L,R) previous to iteration is 3.0858331376860077E-0009 
dist between L and R after applying T 24 is 8. 7084375 703600538E-0005 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 8.7084375703625518E-0005 
iter= 72, dist(L,R) previous to iteration is 8.5043335648104785E-0008 
dist between L and R after applying T 24 is 1.3240893865632102E-0003 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 1.3240893866628222E-0003 
distance between iterates is too large or curvature is big 
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dist between L and R after applying T lb is 2.2719494014277907E-0004 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 2.2719494016150623E-0004 
D dist between iterates continues to be too large or curvature is big 

dist between L and R after applying T 4 is 7.0995114872475991E-0006 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 7.0995114872475992E-0006 

iter= 74, dist(L.R) previous to iteration is 6.9331166869485615E-0009 

dist between L and R after applying T 24 is 1.3422129968050613E-0005 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 1.3422129968054511E-0005 

iter= 624, dist(L,R) previous to iteration is 4.5100485932166933E-0010 
dist between L and R after applying T 10 is 4.4858004094609686E-0007 
dist(L,Ll)+dist(Ll,M)+dist(M,Rl)+dist(Rl,R)= 4.4858004094609686E-0007 
iter= 629, dist(L,R) previous to iteration is 4.3806644611849589E-0010 

Sup distance from fixed point p to point T 1258 (y ) is 0.00022627982 
LI distance from fixed point p to point T 1258 (j/ ) is 0.01670265394 
Euclidean distance from fixed point p to point T 1258 (j/ ) is 0.00141044319 
Euclidean distance from fixed point p to point L is 0.00141021923 
Euclidean distance from fixed point p to point R is 0.00141066714 
Euclidean distance between L and R is 0.00000044858 
angle between Wg(p) and iterated arc LR — 0.00003 radians, angle in degrees is ~ 
angle between vectors (p, L) and (p, R) is 0.00002 angle in degrees is ~ 
angle between vectors (p,T 2 (L)) and (p,T 2 (R)) is 0.00000 radians, angle in degrees is ~ 
angle between vectors (p,T 4 (L)) and (p,T 4 (R)) is 0.00000 radians angle in degrees is ~ 
angle between vectors (p,T 6 (L)) and (p,T 6 (R)) is 0.00000 radians angle in degrees is ~ 
angle between vectors (p,T 8 (L)) and (p,T 8 (R)) is 0.00000 radians angle in degrees is ~ 
angle between vectors (p, T 10 (L)) and (p,T 10 (R)) is 0.00000 radians angle in degrees is ~ 
angle between vectors (p, T 12 (L)) and (p,T 12 (R)) is 0.00001 radians angle in degrees is ~ 
angle between vectors (p,T 14 (L)) and (p,T 14 (R)) is 0.00009 radians angle in degrees is ~ 
Enter gap distance as a real exponent of 2 between 3 and 20 
last exponent used is 11.00000000 
To finish the program enter exponent=0. Chosen exponent = 0.00000000 
Press ENTER to finish the program. 
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